An optimized hybrid Convolutional Perfectly Matched Layer for efficient absorption of electromagnetic waves

Accepted Manuscript An Optimized Hybrid Convolutional Perfectly Matched Layer for Efficient Absorption of Electromagnetic Waves

Amirashkan Darvish, Bijan Zakeri, Nafiseh Radkani

PII: DOI: Reference:

S0021-9991(17)30868-9 https://doi.org/10.1016/j.jcp.2017.11.030 YJCPH 7729

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

15 September 2017 13 November 2017 23 November 2017

Please cite this article in press as: A. Darvish et al., An Optimized Hybrid Convolutional Perfectly Matched Layer for Efficient Absorption of Electromagnetic Waves, J. Comput. Phys. (2017), https://doi.org/10.1016/j.jcp.2017.11.030

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

An Optimized Hybrid Convolutional Perfectly Matched Layer for Efficient Absorption of Electromagnetic Waves Amirashkan Darvish, Bijan Zakeri, Nafiseh Radkani

Abstract— A hybrid technique is studied in order to improve the performance of Convolutional Perfectly Matched Layer (CPML) in the Finite Difference Time Domain (FDTD) medium. This technique combines the first order of Higdon’s annihilation equation as Absorbing Boundary Condition (ABC) with CPML to vanish the Perfect Electric Conductor (PEC) effects at the end of the CPML region. An optimization algorithm is required to find optimum parameters of the proposed absorber. In this investigation, the Particle Swarm Optimization (PSO) is utilized with two separate objective functions in order to minimize the average and peak value of relative error. Using a standard test, the overall performance of the proposed absorber is compared to the original CPML. The results clearly illustrate this method provides approximately 10 dB enhancements in CPML absorption error. The performance, memory and time requirement of the novel absorber, hybrid CPML (H-CPML), was analyzed during 2D and 3D tests and compared to most reported PMLs. The H-CPML requirement of computer resources is similar to CPML and can simply be implemented to truncate FDTD domains. Furthermore, an optimized set of parameters are provided to generalize the hybrid method. Index Terms— Absorbing Boundary Condition, Numerical Analysis, Finite Difference Method, Perfectly Matched Layer, Optimization methods.

I. INTRODUCTION The finite difference time domain (FDTD) method, introduced by Yee [1], is a differential equation solver. It is extensively used in numerical simulation of electromagnetic waves. One of the most significant challenges in the implementation of this solver is finding a technique to terminate its computational domain. Many techniques are introduced to this end, but the efforts can be divided into two distinct areas. The first research area includes analytical Absorbing Boundary Condition (ABC). One of the initial ABCs for acoustic waves, proposed by Engquist and Majda, reveals the acceptable absorption of outgoing waves [2]. Mur implemented the first and second order of Engquist and Majda’s ABC in electromagnetics [3]. By focusing on annihilation operators, Higdon enhanced Engquist and Majda’s ABC for absorption of outgoing waves in specific directions within complex layered media [4, 5]. Ramahi presented the complementary operators method (COM) in order to cancel the boundary reflections [6]. These methods use the tangential fields near boundaries to calculate the equivalent current at the boundaries using Green theorem. The process usually requires storing the time-history tangential fields on the boundary. In addition, most ABCs only operate for particular wave directions. Otherwise stated, they usually trigger an unwanted reflection in other incident directions, especially in near grazing incidences. These drawbacks have made them inefficient for a broad range of electromagnetic problems. High-order ABCs can be more effective. However, the formulation of Higdon’s ABC is not used for higher orders because the original formulations of these ABCs usually involve high-order derivatives. In [7] a new high-order sequence of Non-Reflecting Boundary Condition (NRBC) has been proposed based on the reformulation of Higdon’s wave-product conditions for timedependent waves. This scheme introduces special auxiliary variables that facilitates the use of an NRBC of any arbitrarily highorder. Furthermore, it involves no high derivatives beyond second order. In [8], however, another new auxiliary variable formulation of high-order RBC has been developed. Besides the previous works, there are some other studies about high-order ABCs that can be found in [9-11]. It should be mentioned that despite the acceptable performance of high-order ABCs, these methods often involve complex implementations. The second structure to truncate FDTD lattice (the second area) is a physical layer around the domain of the problem. By specifying the material of the layer, it can physically absorb the waves and behave like the walls of anechoic chambers. Perfectly Amirashkan Darvish is with the Department of Electrical and Computer Engineering, Babol Noshirvani University of Technology, Iran, (email: [email protected]). Bijan Zakeri is with the Faculty of Electrical and Computer Engineering, Babol Noshirvani University of Technology, Iran, (email: [email protected]). Nafiseh Radkani is with the Electrical and Computer Engineering Department, Babol Noshirvani University of Technology, Iran, (email: [email protected]). Corresponding Author: Bijan Gatabi Zakeri ([email protected]), Shariatee St. P.O. Box 47135-484, Post Code 47148-71167, Babol, Iran.

2 Matched Layers (PMLs), first introduced by Berenger [12], are one of the most popular absorbing media. They can provide highly effective truncation boundaries for an FDTD lattice with very low reflection error. In addition, it is possible for PMLs to be located closer to the object under test. However, Berenger’s split-field PML suffers from weakly absorption of evanescent waves and late time reflections [13]. Stretched Coordinate PML (SC-PML) [14], Near PML (NPML) [15], Uniaxial PML (UPML) [16], Auxiliary Differential Equation PML (ADE-PML) [17] and the Convolutional PML (CPML) [18] are some of the most popular PMLs developed to improve various aspects of Berenger’s work. CPML, uses the convolution theorem and the parameters of the Complex Frequency Shifted PML (CFS-PML) [19] to provide a more effective tool which is independent of the medium. Higher order PMLs are other types of PMLs, which have recently been of interest [20-22]. The higher-order PML introduced in [23] is based on unsplit-field SC-PML formulations and the ADE method and requires less memory and computational time in comparison to those in [20-22]. In [24] the nth-order PMLs, based on the unsplit-field formulations and the Z-transform methods, has been presented. The second-order PML is the optimal choice for truncating arbitrary FDTD domains in general cases. However, while these high-order PMLs are the target of many recent studies, CPMLs are more likely to be used in recent FDTD simulations. References [25, 26] provide comprehensive information about various PMLs. The performance of some higher-order ABCs and PMLs has been compared for two-dimensional frequency-domain problems governed by the Helmholtz equation in [27]. The comparisons were based on the ABC order and the number of PML cells. Another study has outlined a different scenario to compare these methods for a dispersive one-dimensional problem in Finite Element Method (FEM) [28]. Then, a parametric study has been performed to analyze the influence of ABC and PML parameters on their accuracy through a wide range of frequencies. In fact, the parametric criterion proposed in [28] is mostly based on frequency. Furthermore, PML absorption depends on various factors, which should be considered in each problem, such as frequency and medium characteristics. In all of the mentioned works, there is a lack of a reliable criterion to propose best PML parameters for a wide range of problems. The Combination of PMLs with ABCs are firstly introduced in [29] for FEM. Such combinations [30, 31] make a great capacity of absorption rather than a single PML or ABC. In [30], the Berenger’s PML is theoretically combined with conformal ABCs for the FEM solver. The combined absorber leads to better performance over the Berenger’s PML in single frequency scattering problems. The combination of unsplit-field PML with three different ABCs is introduced in [31]. It showed that the unsplit-field PML with Ramahi’s ABC could provide more accurate results in FDTD simulations. By extension of using the FDTD technique in full-wave electromagnetic problems, achieving a high level of accuracy in simulations needs an improvement in the performance of the FDTD absorbers. In order to improve the absorption quality, this paper provides a different parametrical and numerical approach for the combination of CPMLs and first order of Higdon’s ABC in FDTD simulator. To this regard, a CPML backed by first order of Higdon’s ABC has been used to avoid extra reflections. The performance of the hybrid absorber, entitled hybrid CPML (H-CPML), has been analyzed by comparing the absorption errors. A different set of new parameters is required to be used for the hybrid structure. To this regard, an optimization algorithm is used in order to achieve optimal performance of the absorber in all simulation times. It is shown that an appropriate evolutionary technique can use the maximum potentials of the absorber by finding the H-CPML optimum parameters for the special conditions of the specific problem. To avoid having a problem specific technique, a useful lookup table is provided as a criterion for choosing the necessary optimized parameters based on the spatial cell size of the desired problem. Results show that the new configuration of H-CPML can provide greater absorption in comparison with the original CPML. In addition, choosing the firstorder of Higdon’s ABC has the advantage of simple implementation and much less usage of computer resources. Therefore, more performing high-order ABCs, such as those proposed in [7] and [8], have not been considered in this paper. Section II provides some preliminaries of Higdon’s annihilation equation and CPML which are required to construct the hybrid layer. In section III, by focusing on the main features of selected absorbers, the hybrid technique is introduced. Then, in section IV, a standard test is outlined to analyze the overall performance of the H-CPML. In section V, the optimization procedure to find the optimum parameters will be described. Section VI proposes a general criterion for choosing parameters. In section VII, H-CPML is applied to a 2D rigorous test in order to analyze its ability in absorption of evanescent waves. Lastly, section VIII presents the H-CPML performance in a 3D realistic medium. All of the results are compared to some original PMLs, such as SC-PML, UPML and CPML. For the mentioned absorbers, the requirement of computer resources is compared in each test. II. REQUIRED PRELIMINARIES This section introduces two absorbers, which are intended to be combined. First, Higdon’s ABC and then, CPML will be described theoretically. By focusing on the important features in the implementation of each absorber, the required basics of the hybrid technique is provided. A. Higdon’s ABC The ABC utilized in this paper is Higdon’s ABC, introduced in [5]. According to the procedure in [32], assume the incident field to the boundaries, propagating in the x − y plane is in the form of a spectrum of plane waves as

3

& & W (r , t ) = ¦ f (ct − kˆi ⋅ r )

(1)

& where c is the speed of light, r is the position vector, and kˆ is the propagating direction vector defined as kˆ = cosα i xˆ + sin α i yˆ

(2)

where α i is the incident angle to the boundary which is unknown in most problems with the range of [ 0, π / 2] . Higdon proposed the following operator to cancel boundary reflections in some desired directions.

§ ∂ c ∂ ·¸ Η (φ j ; j = 1,  , N ) = ∏ N ¨ + j = 1 ¨ ∂ t cos φ ∂ x ¸ j © ¹

(3)

where φ j is a desired direction of the perfect absorption. When (3) applies to the incident waves of (1), for φ j = α i

ª §∂ & c ∂ ·¸º & »W (r , t ) = 0 HW(r , t ) = «∏ N ¨ + « j =1¨© ∂t cosφ j ∂x ¸¹» ¬ ¼

(4)

In practical problems, α i could have a wide range of values. Therefore, choosing optimal values of φ j requires knowledge about the incident behavior. The parameter N determines the order of the operator. N = 1 provides the first-order of Higdon’s ABC with one desirable angle and N = 2 leads to the second-order formulation, which allows selecting two desired angles. By writing the discretized form of (4), it can be observed that the Higdon’s first-order requires storing one-step of time-history fields near the boundaries. Similarly, the second-order requires two steps [32]. Although higher orders of the Higdon’s ABC can cover a wide range of outgoing wave directions, they usually need to store a huge amount of time-history fields in boundaries, causing inefficiency. Moreover, the discretized formulation of the higher orders is harder to use due to their complexity. B. Convolutional PML The convolutional PML, introduced by Roden and Gedney [18], is known as an efficient absorber in the PMLs family. It has three primary advantages over traditional PMLs: higher attenuation of outgoing waves, absorption of evanescent waves and the independency of the medium. In fact, it can be used in various mediums, such as lossy, dispersive, anisotropic and even nonlinear. In addition, its computational storage is the same as PMLs, which makes it more desirable to use. Based on the SC-PML formulations derived in [14], consider the following frequency domain equation

jωεE x + σE x =

1 ∂ 1 ∂ Hz − Hy s y ∂y s z ∂z

(5)

where E x , H y , and H z are x, y, and z directed fields, respectively. The parameters ω , ε , and σ are angular frequency, permittivity, and conductivity of the medium respectively. Parameter s u is the stretched coordinate matrix proposed in [19] su = k u +

σu , (u = x, y , z ) α u + jωε 0

(6)

where σ u is the conductivity in u direction. α u and σ u are real positive numbers and k u ≥ 1 . Equation (5) can be written in the following time domain convolutional form due to the frequency dependent behavior of s u

ε

∂ ∂ ∂ E x + σE x = s y (t ) ∗ H z − s z (t ) ∗ H y ∂t ∂y ∂z

(7)

where * is convolution operator and s u (t ) is the inverse Fourier transform of 1 / s u (ω ) . Using (6) s u (t ) can be written in the

4 following form: s u (t ) =

ξ u (t ) = −

δ (t ) ku

σu ε 0 k u2

e

−

+ ξ u (t )

(8)

· 1 § σu ¨ +αu ¸t ¸ ε 0 ¨© ku ¹

(9)

u(t )

where the δ (t ) and u(t ) are the impulse function and step function, respectively. Equation (7), then can be written as ε

∂E x 1 ∂H z 1 ∂H y +ψ + σE x = − k y ∂y k z ∂z ∂t

hzy

−ψ

hyz

(10)

where ψ hzy and ψ hyz are given by

ψ hzy = ξ y (t ) ∗ ψ hyz = ξ z (t ) ∗

∂H z ∂y

(11)

∂H y

(12)

∂z

The above auxiliary equations can be simplified by using recursive convolution. If the same procedure is applied to other field components, six sets of equations will be derived which can be discretized to the CPML equations. Choosing parameters is one of the most important factors for minimizing the reflection error of CPMLs and other PMLs family. It has been shown that by gradually increasing the loss factor along with PML depth from zero to a maximum value, reflection error is minimized [13]. Also, in [13, 25] a parametric study is provided to select optimum parameters. Parameters σ u , k u , and α u should be slightly changed along the normal direction to the interface by using the following equations [25] §u· ¸ ©d ¹

m

σ u ( u ) = σ max ¨

(13)

§u· k u ( u ) = 1 + ( k max − 1) ⋅ ¨ ¸ ©d ¹

α u (u ) = α

max

§d −u · ¨ ¸ © d ¹

ma

m

(14) (15)

where m and m a are the scaling orders in the polynomial grading and d is the thickness of PML. Equation (13) creates a multistep-discontinuity of the PML conductivity which increases from 0 at x = 0 to σ max at x = d along with the PML depth and leads to a matched media for better absorption of electromagnetic waves. Similarly, in (14), k increases from 1 at x = 0 , the inner boundary of the PML, to k max at x = d , the perfect electric conductor (PEC) outer boundary [16, 33]. The polynomial scaling of α in (15) comes from the phenomena of late-time (low frequency) reflections. The reason for this phenomena is that the tensor coefficient s u has a pole at f = α u 2π ε 0 which can cause a singularity in the equation; especially at low frequencies where α u is small. An appropriate spatial grading can reduce this singularity in all frequencies. To this end, α should be reduced within the PML so that the low frequency spectrum of the waves propagating through the PML can be appropriately attenuated. In this way, (15) leads to a large value of α at the front PML interface and zero value at the end of it. For more details of the update equations and choosing the optimum parameters, readers can refer to [25, 26, 34]. CPML has shown a much better performance in comparison to the conventional ABCs and other kinds of PMLs. Having minimum reflection errors in CPMLs has transformed them to powerful tools in FDTD simulations. In the following, a parametrical method will be proposed to combine the two types of absorbers.

5

Fig. 1. Configuration of absorbers in a 2D TEz FDTD lattice. (a) A medium with Higdon’s ABC (double line) which physically uses only one layer of the FDTD grid. (b) A medium with CPML, which deals with a multi-layer media backed by a perfect electric conductor (PEC with single bold line). (c) A Hybrid CPML, which has been created by substituting the PEC with Higdon’s ABC in CPML media.

III. HYBRID CPML A. Physical Structure One of the primary differences among the PMLs family and the conventional ABCs is that PMLs possess a physical multilayer medium to absorb outgoing waves, however, the first-order Higdon’s ABC, only requires one physical layer to perform as an artificial boundary condition. Having discretized for the TE z mode, the computational FDTD domain of Fig. 1 (a) has been truncated by Higdon’s ABC. This kind of boundary condition is only applied to the latest tangential fields, which are placed on the boundaries and store some time-history fields. As mentioned before, the absorption quality of these absorbers is not sufficiently reliable in all directions of the incident waves. As a result, some unwanted weakened reflections remain in the medium. Fig. 1 (b) demonstrates part of a computational domain comprised of two layers of CPML. The last nodes in the CPML medium do not participate in FDTD update equations. Therefore, they behave like a PEC wall. As illustrated in Fig. 1 (b) This property causes a complete reflection into CPML for waves reaching the end. The reflected waves from PEC continue to attenuation in the reversed direction. The waves that travel into CPML in non-vertical directions are attenuated more than the waves travel in vertical directions. The reason is that in non-vertical cases, the waves travel longer distance into CPML toward PEC and the same distance in the specular directions. In the vertical incidences, it may cause some reflection errors due to the shorter distance the wave travels in CPML. In other words, the duration in which the wave travels into CPML might not be sufficient to cause significant absorption in vertical incidences (this situation usually occurs when the CPML thickness is low). Indeed, the physical nature of CPMLs, which is an absorber layer backed by the PEC wall, can cause some unwanted reflections. In addition to PEC reflections, in the computational domain, there are other unwanted reflections due to the medium discontinuities or the finite spatial sampling (also called discretization error). For a CPML with a fixed thickness, these secondary reflections can be reduced by selecting lower values of σ max with an appropriate scaling order to create small steps of conductivity rising along with PML depth. This approach can cause some of the remaining PEC reflections due to the less attenuation of the outgoing waves, especially in vertical incidences. As a result, PEC reflections can dominate. Therefore, there should be a trade-off between PEC reflections and discretization reflections in CPMLs, which necessitates a difficult decision for setting CPML parameters. By eliminating one of the two causes for reflection error, the other one can be controlled without accepting any trade-off. Fig. 1 (c) shows the third possible configuration for an absorber which consists of a CPML backed by a first order of Higdon’s ABC. The first-order of Higdon’s ABC has been chosen because of its simplicity and less memory usage among other conventional ABCs. Consequently, if any amplitude of incident waves reached the end of the CPML, Higdon ABC will absorb it. In other words, σ max can be further decreased to reduce the discretization error without concerning about the relatively strong remaining waves. Fig. 2 pictorially demonstrates this expected procedure.

6

Fig. 2. A simple graphical model of how a Hybrid CPML is expected to reduce the reflection error without considering the effect of grading order. (a) CPML: choosing causes some reflection error due to discretization reflections. (b) Hybrid CPML: By choosing a value of lower than , the discretization error can be reduced, and the remaining amplitude of waves reaching to the end of the H-CPML can be absorbed by the Higdon’s ABC. Black arrows show the discretization reflected waves. Higher reflections are shown with larger arrows.

B. Numerical Implementation Using the configuration of the hybrid absorber, an improvement in the absorption performance is expected due to the reduction in the PEC reflections. To implement the hybrid absorber, entitled Hybrid CPML (H-CPML), the discretized form of the selected ABC needs to be applied to the last tangential nodes of the FDTD’s medium. Therefore, the following discretized form of the first-order annihilation boundary condition should be used at the end of the CPML region to cancel PEC reflections n +1 E tangential

+

1, j , k

n −1 = E tangential

2, j , k

cos(φ ) ⋅ dx − v ⋅ dt n × ( E tangential cos(φ ) ⋅ dx + v ⋅ dt

1, j , k

n − E tangential

2, j , k

)

(16)

where ϕ , dx, v, and dt are the desired absorption angle, cell size along x direction, wave velocity at the implementation node, and time-steps value, respectively. In this paper, the values of ϕ and v have been fixed to 0 and 3e8, respectively. The second step for H-CPML implementation is selecting the suitable conductivity rising as described in Fig. 2. In order to analyze the proposed technique, a standard test procedure is provided in the next section. The test illuminates that how an evolutionary technique can help an absorber perform with its maximum potential and obtain a new criterion for opting CPML parameters. Then, a comparison will be presented to show the effectiveness of H-CPML.

IV. STANDARD TEST PROCEDURE There are many test procedures, which can effectively analyze the absorption performance of an absorber in FDTD medium. For the first test, consider a problem of finding an absorber’s relative error in a 40 mm × 40 mm TE z FDTD medium as illustrated in Fig. 3. The medium is covered with the desired absorber. A y-directed differentiated Gaussian electric current source is located in the center and given by the following equation

ª t − t0 º J y ( x 0 , y 0 , t ) = −2 × « » ¬ tw ¼

§ t −t · − ¨¨ 0 ¸¸ t e © w ¹

2

(17)

where tw = 26.53 ps and t0 = 4 t w . FDTD space is discretized with Δx = Δy = 1mm , and the Courant factor is set to 0.4950. Electric field ( E y ) is probed in two different points, A and B, and a number of 1200 iterations is considered to achieve the steady states respond. The relative error in each point is defined by the following formula

7

Fig. 3. A 40 mm × 40 mm TEz FDTD medium with two electric field probe at points A and B. The electric current source in center directed toward yaxis. Around the workspace is covered by CPML with thickness d.

ERROR (t ) =

E y (t ) − E yref (t )

(18)

E yref max

ref

where E y (t ) and Ey (t ) are the probed electric fields of the desired points in the 40 mm × 40 mm medium and the electric field ref max

in a 1040 mm × 1040 mm medium with the same distance relative to the central source, respectively. The value of Ey

is the

maximum amplitude of the reference field sensed at the local points. The same problem is studied in [25] where both 6-cell and 10-cell CPML are examined [25, Fig. 7.3b]. In this test, kmax = 1 and αmax = 0.2 are considered for CPML as mentioned in the literature. Third order polynomial grading is used for CPML parameters, and σ max = σ opt is selected to minimize both of the PEC’s reflections and discretization reflections. The parameter, σ opt can be derived from the following formula

σ opt =

0.8(m + 1)

η 0 Δ ε r ,eff μ r ,eff

(19)

where m, η 0 and Δ are the scale order in polynomial grading, free space wave impedance and the lattice cell size, respectively. ε r ,eff and μ r , eff are the effective permittivity and effective permeability of the medium. For this test, both of them are assumed to be 1. Here, 6- and 10-cell of CPML and H-CPML are implemented for analysis. First, the same parameters are considered for both absorbers to see the effect of adding Higdon’s ABC without any change to the CPML construction parameters. Then, considering the H-CPML concept of Fig. 2, lower values of σ max will be tried in order to reduce the discretization error.

V. SIMULATION RESULTS & DISCUSSIONS The first part of this section illuminates the potential of H-CPML in reduction of discretization error, and the second part elaborates the utilization of evolutionary algorithms for opting best parameters. A. Absorption Performance Fig. 4 shows the relative error obtained from (18) at two local points, A and B, for 6-cell and 10-cell of CPML and H-CPML without any changes in construction parameters. There is no difference between the absorption performance of the CPML and HCPML in both cases of 6-cell and 10-cell. Indeed, by choosing the third order of polynomial grading (m = 3) and σ max = σ opt it is guaranteed that significantly small amplitudes of incident waves reach the end of the absorbers. As a result, Higdon’s ABCs may not play an effective role in the absorption. Therefore, it can be inferred that all of the current errors are the result of discretization reflections, and there is no reflection (or probably a significantly small reflection) from PEC to the outside of

8

Fig. 4. The 6-cell and 10-cell relative error of points A and B for CPML and Hybrid CPML with first-order Higdon’s ABC. σ max = σ opt is considered for both types of absorber.

Fig. 5. Relative error of a 6-cell CPML in comparison with the 6-cell HCPML when the value of σmax decreases (Point B).

CPML. Based on this fact and the idea of Fig. 2, it is expected that by choosing smaller values of σ max (σ max < σ opt ) with an appropriate order of grading (m) , any amplitude of a weakened wave that reaches the end of the CPML will be absorbed by Higdon’s ABC. As a result, a decrease in relative error is expectable. A test has been provided in Fig. 5 to analyze this procedure. For the case of 6-cell CPML with the third order of polynomial grading, the value of σ opt is 8.4881. As shown in Fig. 5, for the case of H-CPML construction, by selecting arbitrary values of

σ max (lower than the optimum value), the overall relative error is notably reduced with a slight rise at the beginning of the simulation time. It is obvious that by choosing a very small value for σ max , a weaker absorption occurs into CPML part of the HCPML. As a result, intense waves will reach the end of the H-CPML, and Higdon’s ABC cannot effectively absorb them. This simple test shows the H-CPML capability of absorbing electromagnetic waves. It also implies that the effectiveness of H-CPML is highly dependent on the conductivity distribution. On the other hand, some parameters such as grading order (m) and maximum conductivity ( σ max ) should be selected for utilizing the maximum capability of the absorber. Next subsection shows that an optimization technique helps enhance the H-CPML average or maximum error. B. Optimization procedure As mentioned before, in order to decrease the error level of the whole period and making H-CPML more efficient, an appropriate choice for m and σ max is required. This choice can also help the H-CPML use all of the potentials provided by the Higdon’s ABC. This paper offers PSO algorithm in order to find optimum values of σ max and m for H-CPML [35]. The algorithm can provide a reliable convergence to the optimum performance. Note that by using both of the m and σ max as PSO variables, every value for σ

max

may be possible (less or greater than the σ opt ). To achieve the goal, the following objective

functions can be used: Objective Function 1 = 1 T

³

T 0

2

2

( ERRORA (t) + ERRORB (t) ) dt

(20)

Objective Function 2 = max { ERROR AdB ( t )} + max { ERROR BdB ( t )}

(21)

where in (20), T is the whole period of simulation time and ERROR (t ) is the error function over time, calculated using (18). Objective function of (20) necessitates a decrease in relative error during the whole simulation time, which leads to an improvement in the average error. On the other side, the objective function of (21) provides a relative error with lower maximum

9

ܵ‫ݏݐ݈ݑݏܴ݁ ݁ݎ݋ݐ‬

ܵ݁‫ݏݎ݁ݐ݁݉ܽݎܱܽܲܵܲݐ‬

YES ‫ ܪ݁ݖ݈݅ܽ݅ݐ݅݊ܫ‬െ ‫ܮܯܲܥ‬ ܲܽ‫ݏݎ݁ݐ݁݉ܽݎ‬ ሺ݉ǡ ߪ௠௔௫ ሻ

ܴ‫ݎ݁ݒ݈݋ܵܦܶܦܨ݊ݑ‬ ݅݊‫݉ݑ݅݀݁݉ݐݏ݁ݐ݄݁ݐ‬

No

‫݊݋݅ݐ݅݀݊݋ܥ݀݊ܧ‬

ܷ‫݁ݐܽ݀݌‬ ‫ ܪ‬െ ‫ݏݎ݁ݐ݁݉ܽݎܽܲܮܯܲܥ‬ ܷ‫ݏ݁݅ݐ݅ܿ݋݈ܸ݁݁ݐܽ݀݌‬

‫ݏݏ݁݊ݐݐ݅ܨ݁ݐܽݑ݈ܽݒܧ‬ ܷ‫݃݊݅ݏ‬ ܱܾ݆݁ܿ‫ݏ݊݋݅ݐܿ݊ݑ݂݁ݒ݅ݐ‬

‫ܲ݁ݐ݈ܽݑ݈ܿܽܥ‬஻௘௦௧ ‫݄݂ܿܽ݁݋‬ ܲܽ‫ܩ݈݀݊ܽ݁ܿ݅ݐݎ‬஻௘௦௧ ‫݊݋݅ݐ݈ܽݑ݌݋ܲ ݂݋‬

Fig. 6. PSO block diagram for optimization of H-CPML parameters.

TABLE I THE PSO PARAMETERS Name

cognitive rate ( c1 )

Social rate ( c2 )

Inertial weight (w)

Maximum Velocity ( Vmax )

Value

2.00

2.00

0.9-0.4

4

value; it can guarantee no high reflection in entire simulation time. Table I shows the PSO parameters used for the optimization process. Some previous works [36, 37] have shown that the parameters of Table I are the ideal sets for better convergence of the algorithm in various applications. For PSO, the number of iterations and population size are set to 100 and 10, respectively. In Fig. 6, a flowchart has been provided which describes the optimization procedure for the H-CPML in a specified FDTD problem. Particle locations are substituted by PML parameters, and can be updated over iterations until the stop condition is reached. Fig. 7(a) and (b) illustrate the relative error of points A and B, respectively when the first objective function of (20) is used to optimize the construction parameters of H-CPML (m and σ max ). The PSO algorithm has found m = 2.0729 and

σ max = 5.2808 . By using the optimum parameters, an effective improvement in H-CPML absorption performance is achieved during the entire period, as shown in Fig. 7. It provides a lower level of error in most periods, especially in steady states. As mentioned, (20) causes a suppression in the average error. Therefore, the optimized results of H-CPML in Fig. 7 only have an average improvement in the whole simulation time, not in maximum values. Fig. 8(a) and (b), show the relative error in points A and B, respectively, when the second objective function of (21) has been used as objective function. H-CPML provides -80.5 dB and -79.7 dB maximum error of the points A and B, respectively. While, standard CPML provides -73.5 dB and -72.1 dB for the local points, A and B, respectively. Therefore, for the 6-cell H-CPML, the results show around 7 dB improvement over the conventional CPML. In order to analyze the performance of H-CPML with other number of cells, the same test procedure has been outlined for the absorber with various cell numbers. For each test, 2000 time steps have been considered. PSO is used for each test and the results are depicted in Fig. 9. Fig. 9(a) shows the average values of relative error with respect to the number of cells, when the objective function of (20) has been used for optimizing H-CPML parameters (m and σ max ). For all cell numbers, H-CPML has provided more than 10 dB improvement over original CPML of [25]. Fig. 9(b) also shows the maximum values of relative error versus cell number when the objective function of (21) has been used. Again, H-CPML performs better than the original CPML. In addition, it can be inferred that by increasing the number of cells (N=10 for example) the improvements over the original CPML increases (this increment starts from five cells). The example ensures that H-CPML is capable of showing better performance in general problems, especially when the CPML part consists of 9 or 10 cells. It also provides useful information about the effects of using optimum parameters in H-CPML and its absorption potentials. The following section introduces a criterion to generalize the optimization approach into all problems.

10

Fig. 7. Relative error of a 6-cell H-CPML optimized by PSO algorithm with the first objective function of (20) in compare with the 6-cell of original CPML. (a) Point A (b) Point B

Fig. 8. Relative error of a 6-cell H-CPML optimized by PSO algorithm with the second objective function of (21) in compare with the 6-cell of original CPML. (a) Point A (b) Point B

Fig. 9. The average and maximum value of the relative error (in dB) versus various cell numbers for point A. (a) Average value of the relative error when the first objective function of (20) is used (b) peak value of the relative error when the second objective function of (21) is used.

11 TABLE II H-CPML OPTIMUM PARAMETERS Maximum Optimization Spatial Cell Size m 0.01mm

σ

Average Optimization max

m

σ

max

2.2649

680.6959

2.9033

695.5708

0.1mm

2.2423

69.8992

2.5052

66.1108

1mm

2.2565

7.9606

2.2412

8.2416

10mm

2.6878

0.93432

2.6814

0.9929

VI. PARAMETRIC LOOKUP TABLE By using an optimization procedure to define H-CPML optimum parameters, all absorption potential provided by H-CPML can be utilized effectively. This approach becomes impractical when we deal with large realistic mediums. In addition, since E yref max is unknown in most scattering problems, an independent criterion is required for choosing H-CPML optimum parameters rather than repeating same optimizations for each problem. At this stage, this paper introduces a technique that necessitates few independent optimization processes in order to achieve a global criterion. To this regard, considering the approximated equation of (19), it can be inferred that the optimum conductivity is mainly dependent on medium characteristics ( ε r ,eff and μ r , eff ) and mesh size ( Δ ). As a result, for a specified mesh size and medium characteristic, the optimum values of the parameters (m and σ max ) can be defined separately and used for any other problems with the same spatial cell size and same medium. Assuming the CPML region is in free space (scatterers are usually placed in the middle of simulation medium, out of CPML), the values of ε r ,eff and μ r , eff can be set to 1 in all simulations. In other words, with only one procedure of optimization for a small medium, optimum parameters can be utilized for every other simulation with the same spatial cell size. To obtain the goal, optimum parameters should be extracted from hardest situation in order to be valid in all other problems. Considering this, an arbitrary condensed medium with a randomly shaped scatterer close to the H-CPML boundaries has been considered for various cell sizes. A number of 50 optimization process have been performed for each specified cell size (with its own random scatterer) to reduce the error level in the whole medium. An average of the associated optimum parameters for 10-cell H-CPML are depicted in Table II. To remove the effect of frequency on H-CPML absorption, a wideband Gaussian pulse ranging from zero to maximum valid frequency is considered as the excitation source. As a result, the provided parameters are independent of frequency and problem type. In order to further investigation of the H-CPML performance with the proposed criterion, in the following section, a rigorous test procedure will be outlined. In addition, the H-CPML usage of computer resources, simulation time and absorption performance of evanescent waves will be analyzed. The term ‘rigorous’ is selected because of the specific type of the problem and the intensive conditions provided by decreasing the spatial size of the problem.

VII. RIGOROUS TEST PROCEDURE A. Define Rigorous Problem In the second test, a PEC layer having a thick of 1-Cell is located inside a 110 cell × 7 cell of a 2D FDTD medium. As illustrated in Fig. 10, the absorber walls are too close to the PEC layer (3 cells away from PEC in x and y directions). The PEC length is boundless along z-direction and the width is 104 cells. As more details, the PEC layer is created by considering an extremely large number to both of σ x and σ y in a fixed line ( σ is the medium conductivity). A y-directed current source with a differentiated Gaussian pulse (as stated in the previous test) is located at the center. The FDTD lattice is discretized with Δx = Δy = 1mm and Δt = 1.1785 ps . The -component of electric field is probed at point C, and the relative error is calculated using (18). A number of 1500 iterations is considered to achieve steady states. The interaction of the current source and the PEC layer produces an intensive TE z polarized evanescent wave in the point C. Therefore, an absorber with greater reduction of evanescent waves is preferred for such problems. SC-PML [14], CPML and HCPML are applied to this test. For standard SC-PML, tensor coefficients are assumed as m = 4 , σ max = 0.7 σ opt , kmax = 11 and for CPML, m = 4 , σ max = 0.9 σ opt , k max = 9 , and αmax = 0.05 are chosen as standard parameters. The parameter, σ opt , can be derived from (19). In H-CPML, all of the construction parameters are identical to the CPML, except that m and σ max which are selected using Table II. This makes it possible for Higdon’s ABC to effectively take part in absorption.

12

Fig. 10. Geometrical view of a 110 mm × 7 mm TE z FDTD medium with a PEC layer inside it. The electric current source in center directed toward y-axis. Around the workspace is covered by desired absorber. Y-directed electric field is probed at point C.

Fig. 11. Relative error versus time for 10-cell H-CPML optimized by PSO algorithm using the second objective function of (21) in compare with the 10cell of original CPML, conventional SC-PML and CPML with same optimized parameters.

B. Simulation Results The relative error of the tested absorbers is illustrated in Fig. 11. H-CPML has provided lower relative error in all simulation times. The maximum errors of the SC-PML, CPML and H-CPML are -56.47 dB, -69.64 dB and -79.92 dB, respectively. Therefore, H-CPML has obtained more than 10 dB error enhancement relative to other standard absorbers with their own proposed parameters. In order to compare the required memory and time, the test medium is extended to 1100 cells × 70 cells and the number of iterations has increased to 5000. The required time and memory is depicted in Table III. The required time has averaged over 50 independent runs. The H-CPML time and memory usage is slightly greater than CPML. This is due to the fact that the Higdon’s annihilation equation, which is mounted on the last nodes of the computational domain, only requires to store first step of timehistory fields. Otherwise stated, the first-order of Higdon ABC applies a tiny computational process on system, which results the same computational time for H-CPML. It should be mentioned that the whole test is implemented using 8 processors of Intel core i7 3612QM 2.10 GHz and 8 GB RAM.

VIII. 3D TEST PROCEDURE A. Define 3D Problem In this section, another problem is outlined to show the effectiveness of H-CPML in 3D realistic problems. Consider an elongated thin PEC plate, which is located in the middle of a 3D free space. The dimensions of the PEC sheet are 100 cells × 25 cells × 1 cells. As demonstrated in Fig. 12, the distance between the absorber walls to the PEC is selected to be different in x, y and z directions (5-cell in direction x, 6-cell in direction y and 3-cell in direction z). The spatial medium is discretized using cubic cell spanning Δx = Δy = Δz = 1mm and time-step length is selected to be Δt = 1.9066 ps . The excitation current is similar to the previous examples, except that the parameter, tw , is set to 53 ps. It is located 1mm above one of the corners of the rectangular scatterer. The z-directed electric field is probed at the opposite corner of the sheet where an edge singularity is presented. The

13

Fig. 12. An illustration of the 3D problem. The problem consists of a medium with a 100 mm × 25 mm × 1 mm PEC sheet at the center.

Fig. 13. Relative error versus time for 10-cell H-CPML optimized by PSO algorithm using the objective function of (21) in compare with the 10-cell of original CPML and UPML.

singularity triggers evanescent waves and also creeping waves supported by the plate which require a long-time interaction with the absorber [18, 24, 25]. Therefore, a number of 2000 time steps has been considered for the test. In order to obtain the reference data, the space between the absorber walls and PEC scatterer has been extended 100-cells in all directions. Therefore, a significantly large domain, 225 mm × 300 mm × 200 mm, is achieved as the reference medium. H-CPML, UPML and CPML are applied to this problem in order to compare their absorption performance. For standard UPML, m = 3 , σ max = 0.75 σ opt , and kmax = 15 are selected. For standard CPML, m = 3 , σ max = σ opt , kmax = 15 , and

αmax = 0.24 are assumed to have the best performance as proposed in [25]. H-CPML uses the same parameters of CPML except for m and σ max , which are selected from Table II. Such 3D tests have been used in some other studies [18, 24] in order to compare the performance of various PMLs. B. Simulation Results Fig. 13 shows the relative error of H-CPML in comparison to CPML and UPML (with the original form of choosing parameters as in [25]). Again, H-CPML has provided lower level of reflection error (equal to -90 dB) as opposed to the conventional CPML and UPML (with -80 and -57.5 dB of relative error, respectively). Table IV shows the memory and time requirement for the absorbers under test. The reference medium with the dimensions of 225 mm × 300 mm × 200 mm and the number of 2000 iterations are considered for the test. The same computer system of the previous test was used for simulations. H-CPML running procedure includes extra calculation for the first order of Higdon’s ABC, but the calculations do not significantly affect its performance (as reported in Table IV). Therefore, H-CPML requirements of system resources are similar to conventional CPML. This 3D test ensures that H-CPML can be simply utilized instead of the other common absorbers.

IX. FURTHER DISCUSSIONS In all optimization procedures used in this paper, the parameters k max and α max have been fixed without any change to

14 TABLE III MEMORY AND TIME REQUIREMENT IN THE 2D RIGOROUS TEST Memory Usage (Byte)

Simulation Time (Second)

SC-PML (10 cells)

22695.322K

52.98

CPML (10 cells)

22708.224K

53.09

H-CPML (10 cells)

22827.776K

53.32

SC-PML (16 cells)

30951.652K

61.23

CPML (16 cells)

31086.944K

61.94

H-CPML (16 cells)

31279.168K

62.23

TABLE IV MEMORY AND TIME REQUIREMENT IN THE 3D TEST Memory Usage (Byte)

Simulation Time (Second)

UPML (10 cells)

5010.452M

24751

CPML (10 cells)

5188.188M

25752

H-CPML (10 cells)

5191.135M

25790

UPML (16 cells)

6231.903M

33450

CPML (16 cells)

6949.344M

34619

H-CPML (16 cells)

6955.007M

34682

observe how a good choice of m and σ max can lead to absorption of the remaining waves by Higdon’s ABC, as modeled in Fig 2. Indeed, optimization of the other parameters of H-CPML may result in an improved performance for it, which can be a good topic for future work. In this paper, only the effectiveness of the first order of Higdon’s ABC has been studied. In other words, the combination of CPMLs with other types of the conventional ABCs or the higher orders is another attractive subject for further studies. Although, it should be noted that using higher orders of conventional ABCs requires more computational resources as mentioned before. As we had in CPMLs, H-CPML performance is highly dependent on the excitation source frequency range, the shapes and types of scatterers and medium characteristics [28]. A reliable and ideal criterion for choosing optimum parameters is a criterion that considers all of the mentioned factors. However, by looking at the recent works in this area, it can be observed that there is a lack of a technique, which can provide the best parameters for each specific problem. Using evolutionary algorithms in minimizing the absorber error is a new technique, which can effectively use all of the absorber potentials by including the mentioned factors. However, this technique, in and of itself, suffers from inapplicability in all problems, the H-CPML parametric lookup table provided in this paper is a suitable criterion, which can suggest the best parameters for a wide range of problems.

X. CONCLUSION This paper introduces a new different parametrical and numerical approach for combining CPMLs and first order of Higdon’s annihilation equation in FDTD simulator. The construction consists of a CPML backed by a Higdon’s ABC (entitled H-CPML), which can lead to greater absorption of the remaining waves at the end of CPML. After a parametrical study, the performance of the absorber has been analyzed during three test procedures. The first test has analyzed the overall absorption performance in terms of relative error at two local points. In order to obtain a reliable selection for H-CPML parameters, the PSO algorithm has been utilized for the specified problem. In addition, the performance of H-CPML with various cell numbers has been studied. Then a useful lookup table has been provided for choosing H-CPML optimum parameters. In the second problem, the ability of absorption of evanescent waves has been studied during a 2D rigorous problem in which a PEC scatterer placed too close to absorbers. Our third test has focused on a 3D problem. H-CPML performance has been compared to original UPML and CPML. All provided examples ensure that H-CPML can lead to approximately 10 dB improvement in CPML absorption. Furthermore, the memory and time requirements of H-CPML have been compared to CPML, SC-PML and UPML. The first-order of Higdon’s ABC provides no significant difference in time and memory usage between the H-CPML and CPML, which makes H-CPML more desirable to use.

15

ACKNOWLEDGEMENT This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors

REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]

K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media," IEEE Trans. Antennas Propag, vol. 14, no. 3, pp. 302-307, 1966. B. Engquist and A. Majda, "Absorbing boundary conditions for numerical simulation of waves," Proceedings of the National Academy of Sciences, vol. 74, no. 5, pp. 1765-1766, 1977. G. Mur, "Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations," Electromagnetic Compatibility, IEEE Transactions on, vol. 23, no. 4, pp. 377-382, 1981. R. L. Higdon, "Numerical absorbing boundary conditions for the wave equation," Mathematics of computation, vol. 49, no. 179, pp. 65-90, 1987. R. L. Higdon, "Absorbing boundary conditions for difference approximations to the multidimensional wave equation," Mathematics of computation, vol. 47, no. 176, pp. 437-459, 1986. O. M. Ramahi, "The complementary operators method in FDTD simulations," Antennas and Propagation Magazine, IEEE, vol. 39, no. 6, pp. 33-45, 1997. D. Givoli and B. Neta, "High-order non-reflecting boundary scheme for time-dependent waves," Journal of Computational Physics, vol. 186, no. 1, pp. 2446, 2003. T. Hagstrom and T. Warburton, "A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and extensions to first-order systems," Wave motion, vol. 39, no. 4, pp. 327-338, 2004. D. Givoli and B. Neta, "High-order non-reflecting boundary conditions for dispersive waves," Wave Motion, vol. 37, no. 3, pp. 257-271, 2003. D. Givoli, T. Hagstrom, and I. Patlashenko, "Finite element formulation with high-order absorbing boundary conditions for time-dependent waves," Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 29, pp. 3666-3690, 2006. T. Hagstrom, A. Mar-Or, and D. Givoli, "High-order local absorbing conditions for the wave equation: Extensions and improvements," Journal of Computational Physics, vol. 227, no. 6, pp. 3322-3357, 2008. J.-P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," Journal of computational physics, vol. 114, no. 2, pp. 185-200, 1994. J.-P. Bérenger, "Improved PML for the FDTD solution of wave-structure interaction problems," Antennas and Propagation, IEEE Transactions on, vol. 45, no. 3, pp. 466-473, 1997. W. C. Chew and W. H. Weedon, "A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates," Microwave and optical technology letters, vol. 7, no. 13, pp. 599-604, 1994. S. A. Cummer, "A simple, nearly perfectly matched layer for general electromagnetic media," Microwave and Wireless Components Letters, IEEE, vol. 13, no. 3, pp. 128-130, 2003. S. D. Gedney, "An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices," Antennas and Propagation, IEEE Transactions on, vol. 44, no. 12, pp. 1630-1639, 1996. O. Ramadan, "Auxiliary differential equation formulation: an efficient implementation of the perfectly matched layer," IEEE microwave and wireless components letters, vol. 13, no. 2, pp. 69-71, 2003. J. A. Roden and S. D. Gedney, "Convolutional PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media," Microwave and optical technology letters, vol. 27, no. 5, pp. 334-338, 2000. M. Kuzuoglu and R. Mittra, "Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers," Microwave and Guided Wave Letters, IEEE, vol. 6, no. 12, pp. 447-449, 1996. D. Correia and J. M. Jin, "Performance of regular PML, CFS PML, and second order PML for waveguide problems," Microwave and Optical Technology Letters, vol. 48, no. 10, pp. 2121-2126, 2006. M. W. Chevalier and U. S. Inan, "A PML using a convolutional curl operator and a numerical reflection coefficient for general linear media," Antennas and Propagation, IEEE Transactions on, vol. 52, no. 7, pp. 1647-1657, 2004. D. Correia and J.-M. Jin, "On the development of a higher-order PML," Antennas and Propagation, IEEE Transactions on, vol. 53, no. 12, pp. 4157-4163, 2005. N. Feng and J. Li, "Novel and efficient FDTD implementation of higher-order perfectly matched layer based on ADE method," Journal of Computational Physics, vol. 232, no. 1, pp. 318-326, 2013. N. Feng, Y. Yue, C. Zhu, L. Wan, and Q. H. Liu, "Second-order PML: Optimal choice of nth-order PML for truncating FDTD domains," Journal of Computational Physics, vol. 285, pp. 71-83, 2015. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method. Artech House, 2005. J.-P. Bérenger, "Perfectly matched layer (PML) for computational electromagnetics," Synthesis Lectures on Computational Electromagnetics, vol. 2, no. 1, pp. 1-117, 2007. D. Rabinovich, D. Givoli, and E. Bécache, "Comparison of high order absorbing boundary conditions and perfectly matched layers in the frequency domain," International Journal for Numerical Methods in Biomedical Engineering, vol. 26, no. 10, pp. 1351-1369, 2010. G. Lancioni, "Numerical comparison of high-order absorbing boundary conditions and perfectly matched layers for a dispersive one-dimensional medium," Computer Methods in Applied Mechanics and Engineering, vol. 209, pp. 74-86, 2012. J.-M. Jin and W. C. Chew, "Combining PML and ABC for the finite-element analysis of scattering problems," Microwave and Optical Technology Letters, vol. 12, no. 4, pp. 192-197, 1996. Y. Xiao and Y. Lu, "Combination of PML and ABC for scattering problem," IEEE transactions on magnetics, vol. 37, no. 5, pp. 3510-3513, 2001. N. V. Kantartzis and T. D. Tsiboukis, "A higher-order FDTD technique for the implementation of enhanced dispersionless perfectly matched layers combined with efficient absorbing boundary conditions," IEEE transactions on magnetics, vol. 34, no. 5, pp. 2736-2739, 1998. S. D. Gedney, "Introduction to the finite-difference time-domain (FDTD) method for electromagnetics," Synthesis Lectures on Computational Electromagnetics, vol. 6, no. 1, pp. 1-250, 2011. S. D. Gedney, "The perfectly matched layer absorbing medium," Advances in Computational Electrodynamics: The Finite-Difference Time-Domain Method, pp. 263-344, 1998. A. Z. Elsherbeni and V. Demir, The Finite-difference Time-domain Method for Electromagnetics with MATLAB® Simulations. SciTech Pub., 2009. J. Kennedy and R. Eberhart, "Particle swarm optimization," in Neural Networks, 1995. Proceedings., IEEE International Conference on, 1995, vol. 4, pp. 1942-1948 vol.4. R. C. Eberhart and Y. Shi, "Particle swarm optimization: developments, applications and resources," in evolutionary computation, 2001. Proceedings of the 2001 Congress on, 2001, vol. 1, pp. 81-86: IEEE.

16 [37] J. Robinson and Y. Rahmat-Samii, "Particle swarm optimization in electromagnetics," Antennas and Propagation, IEEE Transactions on, vol. 52, no. 2, pp. 397-407, 2004.