# Accurate and Efficient Evaluation of Characteristic Modes

###### Abstract

A new method to improve the accuracy and efficiency of characteristic mode (CM) decomposition for perfectly conducting bodies is presented. The method uses the expansion of the Green dyadic in spherical vector waves. This expansion is utilized in the method of moments (MoM) solution of the electric field integral equation (EFIE) to factorize the real part of the impedance matrix. The factorization is then employed in the computation of characteristic modes, which improves the accuracy as well as the computational speed. An additional benefit is a rapid computation of far fields. The method can easily be integrated into existing method of moments (MoM) solvers. Several structures are investigated illustrating the improved accuracy and performance of the new method.

## I Introduction

The method of moments (MoM) solution to electromagnetic field integral equations was introduced by Harrington [1] and has prevailed as a standard in solving open (radiating) electromagnetic problems [2]. While memory-demanding, MoM represents operators as matrices (notably the impedance matrix [1]) allowing for direct inversion and modal decompositions [3]. The latter option is becoming increasingly popular, mainly due to characteristic mode (CM) decomposition [4], a leading formalism in antenna shape and feeding synthesis [5, 6], determination of optimal currents [7, 8], and performance evaluation [9].

Utilization of CM decomposition is especially efficient when dealing with electrically small antennas [10], particularly if they are made solely of perfect electric conductor (PEC), for which only a small number of modes are needed to describe their radiation behavior. Yet, the real part of the impedance matrix is indefinite as it is computed with finite precision [11, 12].
The aforementioned deficiency is resolved in this paper by a two-step procedure. First, the real part of the impedance matrix is constructed using spherical wave expansion of the dyadic Green function [13]. This makes it possible to decompose the real part of the impedance matrix as a product of a spherical modes projection matrix with its hermitian conjugate. The second step consists of reformulating the modal decomposition so that only the standalone spherical modes projection matrix is involved preserving the numerical dynamics^{1}^{1}1The numerical dynamic is defined as the largest characteristic eigenvalue..

The proposed method significantly accelerates the computation of CMs as well as of the real part of the impedance matrix. Moreover, it is possible to recover CMs using lower precision floating point arithmetic, which reduces memory use and speeds up arithmetic operations if hardware vectorization is exploited [3]. An added benefit is the efficient computation of far field patterns using spherical vector harmonics.

The projection on spherical waves in the proposed method introduces several appealing properties. First is an easy monitoring of the numerical dynamics of the matrix, since the different spherical waves occupy separate rows in the projection matrix. Second is the possibility to compute a positive semidefinite impedance matrix which plays important role in an optimal design [14, 8]. A final benefit is the superposition of modes. [6].

The paper is organized as follows. The construction of the impedance matrix using classical procedure is briefly reviewed in Section II-A and the proposed procedure is presented in Section II-B. Numerical aspects of evaluating the impedance matrix are discussed in Section II-C. In Section III, the spherical modes projection matrix is utilized to reformulate modal decomposition techniques, namely the evaluation of radiation modes in Section III-A and CMs in Section III-B. These two applications cover both the standard and generalized eigenvalue problems. The advantages of the proposed procedure are demonstrated on a series of practical examples in this section. Various aspects of the proposed method are discussed in Section IV and the paper is concluded in Section V.

## Ii Evaluation of Impedance Matrix

This paper investigates mode decompositions for PEC structures in free space. The time-harmonic quantities under the convention , with being the angular frequency, are used throughout the paper.

### Ii-a Method of Moments Implementation of the Efie

Let us consider the electric field integral equation (EFIE) [1] for PEC bodies, defined as

(1) |

with being the impedance operator, the incident electric field [15], the current density, the imaginary unit, and the unit normal vector to the PEC surface. The EFIE (1) is explicitly written as

(2) |

where , is the wave number, the free space impedance, and the dyadic Green function for the electric field in free-space defined as [13, 16]

(3) |

Here, is the identity dyadic, and , are the source and observation points. The EFIE (2) is solved with the MoM by expanding the current density into real-valued basis functions as

(4) |

and applying Galerkin testing procedure [16, 17]. The impedance operator is expressed as the impedance matrix , where is the resistance matrix, and the reactance matrix. The elements of the impedance matrix are

(5) |

### Ii-B Spherical Wave Expansion of the Green Dyadic

The Green dyadic (3) that is used to compute the impedance matrix can be expanded in spherical vector waves as

(6) |

where and if , and and if . The regular and outgoing spherical vector waves [18, 19, 20, 13] are and , see Appendix B. The mode index for real-valued vector spherical harmonics is [20, 21]

(7) |

with , , , for even azimuth functions (), and for odd azimuth functions (). Inserting the expansion of the Green dyadic (6) into (5), the impedance matrix becomes

(8) |

For a PEC structure the resistive part of (8) can be factorized as

(9) |

where is used. Reactance matrix, , cannot be factorized in a similar way as two separate spherical waves occur.

Resistance matrix can be written in matrix form as

(10) |

where is the matrix transpose. Individual elements of the matrix are

(11) |

and the size of the matrix is , where

(12) |

is the number of spherical modes and the highest order of spherical mode, see Appendix B. For complex-valued vector spherical harmonics [20] the transpose in (10) is replaced with the hermitian transpose . The individual integrals in (8) are in fact related to the T-matrix method [22, 19], where the incident and scattered electric fields are expanded using regular and outgoing spherical vector waves, respectively. The factorization (6) is also used in vector fast multipole algorithm [23].

The radiated far-field can conveniently be computed using spherical vector harmonics

(13) |

where are the spherical vector harmonics, see Appendix B. The expansion coefficients are given by

(14) |

where the column matrix contains the current density coefficients . The total time-averaged radiated power of a lossless antenna can be expressed as a sum of expansion coefficients

(15) |

### Ii-C Numerical Considerations

The spectrum of the matrices and differ considerably [8, 12]. The eigenvalues of the matrix decrease exponentially and the number of eigenvalues are corrupted by numerical noise, while this is not the case for the matrix . As a result, if the matrix is used in an eigenvalue problem, only a few modes can be extracted. This major limitation can be overcome with the use of the matrix in (11), whose elements vary several order of magnitude, as the result of the increased order of spherical modes with increasing row number. If the matrix is directly computed with the matrix product (10) or equivalently from matrix produced by (5) small values are truncated due to floating-point arithmetic^{2}^{2}2As an example to the loss of significance in double precision arithmetic consider the sum . [24, 25]. Subsequently, the spectrum of the matrix should be computed from the matrix as presented in Section III.

The matrix also provides a low-rank approximation of the matrix , which is the result of the rapid convergence of regular spherical waves. In this paper, the number of used modes in (6) is truncated using a modified version of the expression in [26]

(16) |

where is the highest order of spherical mode, is the radius of the sphere enclosing the scatterer, and is the ceiling function. The resulting accuracy in all treated cases is satisfactory. The order of spherical modes can be modified to trade between accuracy and computational efficiency, where increasing improves the accuracy. Fig. 1 shows the convergence of the matrix for Example R2.

Substitution of the spherical vector waves, introduced in Section II-B, separates (5) into two separate surface integrals reducing computational complexity. Table I presents computation times^{3}^{3}3Computations are done on a workstation with i7-3770 CPU @ GHz and GB RAM, operating under Windows 7. of different matrices^{4}^{4}4Computation time for the matrix is omitted as it takes longer than the matrix , due to Green function singularity. , , , and for the examples given in Table II. As expected, the matrix requires the most computational resources, as it includes both the matrix and . The computation of the matrix using MoM is faster than the matrix since the underlying integrals are regular. The computation of the matrix using (10) takes the least amount of time for most of the examples. The computational gain is notable for structures with more degrees-of-freedom (d-o-f), .

## Iii Modal Decomposition With the Matrix

Modal decomposition using the matrix is applied to two structures; a spherical shell of radius , and a rectangular plate of length and width (App. D), are presented in Table II. Both structures are investigated for different number of d-o-f, RWG functions [29] are used as the basis functions . The matrices used in modal decomposition have been computed using in-house solvers AToM [30] and IDA [31], see Appendix A for details. Results from the commercial electromagnetic solver FEKO [32] are also presented for comparison. Computations that require a higher precision than the double precision arithmetic are performed using the mpmath Python library [27], and the Advanpix Matlab toolbox [33].

Structure | Example | |||
---|---|---|---|---|

S1 | ||||

S2 | ||||

Spherical shell | S3 | |||

Fig. 9 | S4 | |||

S5 | ||||

Rectangular plate | R1 | |||

Fig. 10 | R2 | |||

R3 | ||||

R4 | ||||

Helicopter | H1 | |||

H2 |

Number of properly calculated modes | |||||

Example | |||||

(see Table II) | (17) | (19) | (20) | (24) | |

S2 | |||||

S3 | |||||

S5 | |||||

R1 | |||||

R2 | |||||

R4 |

### Iii-a Radiation Modes

The eigenvalues for the radiation modes [35] are easily found using the eigenvalue problem

(17) |

where are the eigenvalues of the matrix , and are the eigencurrents. The indefiniteness of the matrix poses a problem in the eigenvalue decomposition (17) as illustrated in [8, 12]. In this paper we show that the indefiniteness caused by the numerical noise can be bypassed using the matrix . We start with the singular value decomposition (SVD) of the matrix

(18) |

where and are unitary matrices, and is a diagonal matrix containing singular values of matrix . Inserting (10), (18) into (17) and multiplying from the left with yields

(19) |

where the eigenvectors are rewritten as , and the eigenvalues are . A comparison of procedure (17) and (19) is shown in Table III. For high order , the classical procedure (17) with double numerical precision yields in unphysical modes with negative eigenvalues (negative radiated power) or with incorrect current profile (as compared to the use of quadruple precision). Using double precision, the number of modes which resemble physical reality (called “properly calculated modes” in Table III) is much higher^{5}^{5}5Quantitatively, the proper modes in Table III are defined as those having less than deviation in eigenvalue as compared to the computation with quadruple precision. for the new procedure (19). It is also worth mentioning that the new procedure, by design, always gives positive eigenvalues .

### Iii-B Characteristic Modes (CMs)

The generalized eigenvalue problem (GEP) with the matrix on the right hand side, i.e., serving as a weighting operator [36], is much more involved as the problem cannot be completely substituted by the SVD. Yet, the SVD of the matrix in (18) plays an important role in CM decomposition.

The CM decomposition is defined as

(20) |

which is known to suffer from the indefiniteness of the matrix [12], therefore delivering only a limited number of modes. The first step is to represent the solution in a basis of singular vectors by substituting the matrix in (20) as (10), with (18) and multiplying (20) from the left by the matrix

(21) |

Formulation (21) can formally be expressed as a GEP with an already diagonalized right hand side [37]

(22) |

i.e., , , and .

Since the matrix is in general rectangular, it is crucial to take into account cases where , (12). This is equivalent to a situation in which there are limited number of spherical projections to recover the CMs. Consequently, only limited number of singular values exist. In such a case, the procedure similar to the one used in [11] should be undertaken by partitioning (22) into two linear systems

(23) |

where , , and . The Schur complement is obtained by substituting the second row of (23) into the first row

(24) |

with expansion coefficients of CMs defined as

(25) |

As far as the matrices and in (18) are unitary, the decomposition (22) yields CMs implicitly normalized to

(26) |

which is crucial since the standard normalization cannot be used without decreasing the number of significant digits. In order to demonstrate the use of (24), various examples from Table II are calculated and compared with the conventional approach (20).

The CMs of the spherical shell from Example S2 are calculated and shown as absolute values in logarithmic scale in Fig. 2. It is shown that the number of the CMs calculated by classical procedure (FEKO, AToM) is limited to the lower modes, especially considering the degeneracy of the CMs on the spherical shell [12]. The number of properly found CMs is significantly higher when using (24) than the conventional approach (20) and the numerical dynamic is doubled. Notice that, even (20) where the matrix calculated from (10) yields slightly better results than the conventional procedure. This fact is confirmed in Fig. 3 dealing with Example R2, where the multiprecision package Advanpix is used as a reference. The same calculation illustrates that the matrix contains all information to recover the same number of modes as (24), but this can be done only at the expense of higher computation time^{6}^{6}6For Example S2 the computation time of CMs with quadruple precision is approximately hours..

While (24) preserves the numerical dynamics, the computational efficiency is not improved due to the matrix multiplications to calculate the term in (23). An alternative formulation that improves the computational speed is derived by replacing the matrix with (10) in (20)

(27) |

and multiplying from the left with

(28) |

The formulation (28) is a standard eigenvalue problem and can be written as

(29) |

where , , and . As an intermediary step, the matrix is computed, which is later used to calculate the characteristic eigenvectors . The eigenvalue problem (29) is solved in the basis of spherical vector waves, , that results in a matrix . For problems with the eigenvalue problem is solved rapidly compared with (20) and (24). The computation times for various examples are presented in Table IV for all three formulations where a different number of CMs are compared. For Example H1 the computation time is investigated for the first and modes. The acceleration using (29) is approximately and times when compared with the conventional method (20). The first characteristic mode of Example H1 is illustrated in Fig. 4.

Example | Time to calculate CMs (s) | |||
---|---|---|---|---|

(see Table II) | (20) | (24) | (29) | |

S1 | ||||

S2 | ||||

S4 | ||||

S4 | ||||

R1 | ||||

R3 | ||||

H1 | ||||

H1 | ||||

H2 |

Two tests proposed in [12] are performed to validate the conformity of characteristic current densities and the characteristic far fields with the analytically known values. The results of the former test are depicted in Fig. 5 for Example S2 and S5 that are spherical shells with two different d-o-f. Similarity coefficients are depicted both for the CMs using the matrix (20) and for the CMs calculated by (24). The number of valid modes correlates well with Table III and the same dependence on the quality and size of the mesh grid as in [12] is observed.

Qualitatively the same behavior is also observed in the latter test, depicted in Fig. 6, where similarity of characteristic far fields is expressed by coefficient [12]. These coefficients read

(30) |

where has been evaluated using (14). The results for characteristic far fields computed from the conventional procedure (20) and the procedure presented in this paper (24) are illustrated in Fig. 6.

Lastly, the improved accuracy of using (24) over (20), is demonstrated in the Fig. 7 which shows current profiles, corresponding to a rectangular plate (Example R2), of a collection of the first 30 modes. It can be seen that for modes with high eigenvalues (numerically saturated regions in Fig. 3) the surface current density in left panel, calculated via (20), shows numerical noise, while the evaluation via (24) still yields a correct current profile.

### Iii-C Restriction to TM/TE modes

Matrix , described in Section II-B, contains projections onto TE and TM spherical waves in its odd () and even rows (), respectively. The separation of TE and TM spherical waves can be used to construct resistance matrices and , where only odd and even rows of matrix are used to evaluate (10).

Matrices and can be used in optimization, e.g., in such a case when the antennas have to radiate TM-modes only [38]. With this feature, characteristic modes consisting of only TM (or TE) modes can easily be found. This is shown in Fig. 8, in which the spherical shell (Example S2) and rectangular plate (Example R2) are used to find only TM (capacitive) and TE (inductive) modes, respectively. In case of a spherical shell this separation could have been done during the post-processing. For a generally shaped body this separation however represents a unique feature of the proposed method.

## Iv Discussion

Important aspects of the utilization of the matrix are discussed under the headings implementation aspects, computational aspects and potential improvements.

### Iv-a Implementation Aspects

Unlike the reactance matrix , the resistance matrix suffers from high condition number. Therefore, the combined approach to evaluate the impedance matrix (matrix using matrix , matrix using conventional Green function technique with double integration) takes advantage of both methods and is optimal for, e.g., modal decomposition techniques dealing with the matrix (radiation modes [35], CMs, energy modes [35, 39], and solution of optimization problems [38]). Evaluation and the SVD of the matrix are also used to estimate number of modes, cf. number of modes of the matrix found by (18) and number of CMs found by (24) in Table III.

### Iv-B Computational Aspects

Computational gains of the proposed method are seen in Table I for the matrix and Table IV for the CMs. The formulation (29) significantly accelerates CMs computation when compared with the classical GEP formulation (20). Moreover, it is possible to employ lower precision floating point arithmetic, e.g. float, to compute as many modes as the conventional method that employs higher precision floating point arithmetic, e.g. double. In modern hardware, this can provide additional performance boosts if vectorization is used.

An advantage of the proposed method is that the matrix is rectangular for , allowing independent selection of the parameters and . While the parameter controls the details in the model, the parameter (or alternatively ) controls the convergence of the matrix and the number of modes to be found. In this paper (16) is used to determine the highest spherical wave order for a given electrical size . The parameter can be increased for improved accuracy or decreased for computational gain depending on the requirements of the problem. Notice that the parameter is limited from below by the convergence and the number of desired modes, but also from above since the spherical Bessel function in decays rapidly with as

(31) |

The rapid decay can be observed in Fig. 1, where the convergence of the matrix to double precision for requires only while (16) gives a conservative number of .

### Iv-C Potential Improvements

Even though the numerical dynamic is increased, it is strictly limited and it presents an inevitable, thus fundamental, bottleneck of all modal methods involving radiation properties. The true technical limitation is, in fact, the SVD of the matrix . A possible remedy is the use of high-precision packages that come at the expense of markedly longer computation times and the necessity of performing all subsequent operations in the same package to preserve high numerical precision.

The second potential improvement relies on higher-order basis functions, which can compensate a poor-meshing scheme (that is sometimes unavoidable for complex or electrically large models). It can also reduce the number of basis function so that the evaluation of CMs is further accelerated.

## V Conclusion

Evaluation of the discretized form of the EFIE impedance operator, the impedance matrix, has been reformulated using projection of vector spherical harmonics onto a set of basis functions. The key feature of the proposed method is the fact that the real part of the impedance matrix can be written as a multiplication of the spherical modes projection matrix with itself. This feature accelerates modal decomposition techniques and doubles the achievable numerical dynamics. The results obtained by the method can also be used as a reference for validation and benchmarking.

It has been shown that the method has notable advantages, namely the number of available modes can be estimated prior to the decomposition and the convergence can be controlled via the number of basis functions and the number of projections. The normalization of generalized eigenvalue problems with respect to the product of the spherical modes projection matrix on the right hand side are implicitly done. The presented procedure finds its use in various optimization techniques as well. It allows for example to prescribe the radiation pattern of optimized current by restricting the set of the spherical harmonics used for construction of the matrix.

The method can be straightforwardly implemented into both in-house and commercial solvers, improving thus their performance and providing antenna designers with more accurate and larger sets of modes.

## Appendix A Used Computational Electromagnetics Packages

### A-a Feko

FEKO (ver. 14.0-273612, [32]) has been used with a mesh structure that was imported in NASTRAN file format [40]: CMs and far fields were chosen from the model tree under requests for the FEKO solver. Data from FEKO were acquired using *.out, *.os, *.mat and *.ffe files. The impedance matrices were imported using an in-house wrapper [31]. Double precision was enabled for data storage in solver settings.

### A-B AToM

AToM (pre-product ver., CTU in Prague, [30]) has been used with a mesh grid that was imported in NASTRAN file format [40], and simulation parameters were set to comply with the data in Table II. AToM uses RWG basis functions with the Galerkin procedure [29]. The Gaussian quadrature is implemented according to [34] and singularity treatment is implemented from [41]. Built-in Matlab functions are utilized for matrix inversion and decomposition. Multiprecision package Advanpix [33] is used for comparison purposes.

### A-C Ida

IDA (in-house, Lund University, [31]) has been used with the NASTRAN mesh and processed with the IDA geometry interpreter. IDA solver is a Galerkin type MoM implementation. RWG basis functions are used for the current densities. Numerical integrals are performed using Gaussian quadrature [34] for non-singular terms and the DEMCEM library [42, 43, 44, 45] for singular terms. Intel MKL library [28] is used for linear algebra routines. The matrix computation routines are parallelized using OpenMP [46]. Multiprecision computations were done with the mpmath Python library [27].

## Appendix B Spherical Vector Waves

General expression of the (scalar) spherical modes is [13]

(32) |

with and being the wavenumber. The indices are , and [20, 21]. For regular waves is a spherical Bessel function of order , irregular waves is a spherical Neumann function, and are spherical Hankel functions for the ingoing and outgoing waves, respectively. Spherical harmonics are defined as [13]

(33) |

with the Neumann factor, the Kronecker delta function and the normalized associated Legendre functions [47].

The spherical vector waves are [20, 13]

(34a) | ||||

(34b) |

where are the radial function of order defined as

, | (35a) | ||||

, | (35b) | ||||

, | (35c) |

with and denotes the real-valued vector spherical harmonics defined as

(36a) | ||||

(36b) |

where denotes the ordinary spherical harmonics [13]. The radial functions can be seperated into real and imaginary parts as

(37) | ||||

(38) |

## Appendix C Associated Legendre Polynomials

The associated Legendre functions are defined [48] as

(39) |

with

(40) |

being the associated Legendre polynomials of degree and . One useful limit when computing the vector spherical harmonics is [13]

(41) |

The normalized associated Legendre function , is defined as follows

(42) |

The derivative of the normalized associated Legendre function is required when computing the spherical harmonics, and is given by the following recursion relation

(43) |

where , .

## Appendix D Spherical Shell and Rectangular Plate

## Appendix E Radiation Modes

Eigenvalues of the radiation modes for Example S2 and R2 are presented in Fig. 11 and Fig. 12. The eigenvalues are computed using both the conventional (17) and the proposed (19) method. It can be seen that the number of modes computed using (19) is significantly higher compared to (17) for both examples. Eigenvalues calculated using quadruple precision SVD of the matrix are also included. The number of correct radiation modes is shown in Table III.

If eigenvalues of the different mesh grids are to be compared the MoM matrices must be normalized. The normalized matrices are , , , , where is the diagonal matrix of basis functions’ reciprocal edge lengths, i.e., .

## References

- [1] R. F. Harrington, Field Computation by Moment Methods. Wiley – IEEE Press, 1993.
- [2] M. N. O. Sadiku, Numerical Techniques in Electromagnetics with Matlab, 3rd ed. CRC Press, 2009.
- [3] G. H. Golub and C. F. Van Loan, Matrix Computations. Johns Hopkins University Press, 2012.
- [4] R. F. Harrington and J. R. Mautz, “Theory of characteristic modes for conducting bodies,” IEEE Trans. Antennas Propag., vol. 19, no. 5, pp. 622–628, 1971.
- [5] B. Yang and J. J. Adams, “Systematic shape optimization of symmetric mimo antennas using characteristic modes,” IEEE Trans. Antennas Propag., vol. 64, no. 7, pp. 2668–2678, July 2016.
- [6] M. Capek, P. Hazdra, and J. Eichler, “A method for the evaluation of radiation Q based on modal approach,” IEEE Trans. Antennas Propag., vol. 60, no. 10, pp. 4556–4567, Oct. 2012.
- [7] M. Capek and L. Jelinek, “Optimal composition of modal currents for minimal quality factor Q,” IEEE Trans. Antennas Propag., vol. 64, no. 12, pp. 5230–5242, 2016.
- [8] M. Gustafsson, D. Tayli, C. Ehrenborg, M. Cismasu, and S. Norbedo, “Antenna current optimization using MATLAB and CVX,” FERMAT, vol. 15, no. 5, pp. 1–29, May–June 2016. [Online]. Available: http://www.e-fermat.org/articles/gustafsson-art-2016-vol15-may-jun-005/
- [9] M. Vogel, G. Gampala, D. Ludick, U. Jakobus, and C. Reddy, “Characteristic mode analysis: Putting physics back into simulation,” IEEE Antennas Propag. Mag., vol. 57, no. 2, pp. 307–317, April 2015.
- [10] Y. Chen and C.-F. Wang, “Electrically small UAV antenna design using characteristic modes,” IEEE Trans. Antennas Propag., vol. 62, no. 2, pp. 535–545, Feb. 2014.
- [11] R. F. Harrington and J. R. Mautz, “Computation of characteristic modes for conducting bodies,” IEEE Trans. Antennas Propag., vol. 19, no. 5, pp. 629–639, Sept. 1971.
- [12] M. Capek, V. Losenicky, L. Jelinek, and M. Gustafsson, “Validating the characteristic modes solvers,” IEEE Trans. Antennas Propag., vol. 65, no. 8, pp. 4134–4145, 2017.
- [13] G. Kristensson, Scattering of Electromagnetic Waves by Obstacles. Edison, NJ: SciTech Publishing, an imprint of the IET, 2016.
- [14] M. Gustafsson and S. Nordebo, “Optimal antenna currents for Q, superdirectivity, and radiation patterns using convex optimization,” IEEE Trans. Antennas Propag., vol. 61, no. 3, pp. 1109–1118, 2013.
- [15] R. F. Harrington, Time-Harmonic Electromagnetic Fields, 2nd ed. Wiley – IEEE Press, 2001.
- [16] W. C. Chew, M. S. Tong, and B. Hu, Integral Equation Methods for Electromagnetic and Elastic Waves. Morgan & Claypool, 2009.
- [17] W. C. Gibson, The Method of Moments in Electromagnetics. CRC press, 2014.
- [18] P. M. Morse and H. Feshbach, Methods of Theoretical Physics. New York, NY: McGraw-Hill, 1953, vol. 2.
- [19] P. C. Waterman, “Symmetry, unitarity, and geometry in electromagnetic scattering,” Phys. Rev. D, vol. 3, no. 4, pp. 825–839, 1971.
- [20] J. E. Hansen, Ed., Spherical Near-Field Antenna Measurements, ser. IEE electromagnetic waves series. Stevenage, UK: Peter Peregrinus Ltd., 1988, no. 26.
- [21] M. Gustafsson and S. Nordebo, “Characterization of MIMO antennas using spherical vector waves,” IEEE Trans. Antennas Propag., vol. 54, no. 9, pp. 2679–2682, 2006.
- [22] P. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE, vol. 53, no. 8, pp. 805–812, Aug. 1965.
- [23] Y. G. Liu, W. C. Chew, L. Jiang, and Z. Qian, “A memory saving fast A-EFIE solver for modeling low-frequency large-scale problems,” Applied Numerical Mathematics, vol. 62, no. 6, pp. 682–698, 2012.
- [24] D. Zuras, M. Cowlishaw, A. Aiken, M. Applegate, D. Bailey, S. Bass, D. Bhandarkar, M. Bhat, D. Bindel, S. Boldo et al., “IEEE standard for floating-point arithmetic,” IEEE Std 754-2008, pp. 1–70, 2008.
- [25] R. Burden, J. Faires, and A. Burden, Numerical Analysis. Cengage Learning, 2015.
- [26] J. Song and W. C. Chew, “Error analysis for the truncation of multipole expansion of vector green’s functions [em scattering],” IEEE Microwave and Wireless Components Letters, vol. 11, no. 7, pp. 311–313, 7 2001.
- [27] F. Johansson et al. (2013, December) mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.18). [Online]. Available: http://mpmath.org/
- [28] Intel. (2017) Intel Math Kernel Library 2017 update 3. [Online]. Available: https://software.intel.com/en-us/mkl
- [29] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propag., vol. 30, no. 3, pp. 409–418, May 1982.
- [30] (2017) Antenna Toolbox for MATLAB (AToM). Czech Technical University in Prague. [Online]. Available: www.antennatoolbox.com
- [31] D. Tayli. (2017) IDA (Integrated Development toolset for Antennas). Lund University.
- [32] Altair. (2016) FEKO. Altair. [Online]. Available: www.feko.info
- [33] Advanpix. (2016) Multiprecision Computing Toolbox for MATLAB. [Online]. Available: http://www.advanpix.com/
- [34] D. A. Dunavant, “High degree efficient symmetrical gaussian quadrature rules for the triangl,” International Journal for Numerical Methods in Engineering, vol. 21, pp. 1129–1148, 1985.
- [35] K. R. Schab and J. T. Bernhard, “Radiation and energy storage current modes on conducting structures,” IEEE Trans. Antennas Propag., vol. 63, no. 12, pp. 5601–5611, Dec. 2015.
- [36] J. H. Wilkinson, The Algebraic Eigenvalue Problem. Oxford University Press, 1988.
- [37] G. Angiulli and F. Venneri, “Use of the simultaneous diagonalization technique in the eigenproblem applied to the computation of the characteristic modes,” ACES Journal, vol. 17, no. 3, pp. 232–238, Nov. 2002.
- [38] M. Capek, M. Gustafsson, and K. Schab, “Minimization of antenna quality factor,” IEEE Trans. Antennas Propag., vol. 65, no. 8, pp. 4115â–4123, 2017.
- [39] L. Jelinek and M. Capek, “Optimal currents on arbitrarily shaped surfaces,” IEEE Trans. Antennas Propag., vol. 65, no. 1, pp. 329–341, Jan. 2017.
- [40] (2017) MSC NASTRAN. [Online]. Available: http://www.mscsoftware.com/support/
- [41] T. F. Eibert and V. Hansen, “On the calculation of potential integrals for linear source distributions on triangular domains,” IEEE Trans. Antennas Propag., vol. 43, no. 12, pp. 1499–1502, Dec. 1995.
- [42] A. G. Polimeridis. (2010) Direct evaluation method in computational electromagnetics (DEMCEM). [Online]. Available: https://github.com/thanospol/DEMCEM
- [43] A. G. Polimeridis and T. V. Yioultsis, “On the direct evaluation of weakly singular integrals in galerkin mixed potential integral equation formulations,” IEEE Trans. Antennas Propag., vol. 56, no. 9, pp. 3011–3019, 2008.
- [44] A. G. Polimeridis and J. R. Mosig, “Complete semi-analytical treatment of weakly singular integrals on planar triangles via the direct evaluation method,” International journal for numerical methods in engineering, vol. 83, no. 12, pp. 1625–1650, 2010.
- [45] ——, “On the direct evaluation of surface integral equation impedance matrix elements involving point singularities,” IEEE Trans. Antennas Propag., vol. 10, pp. 599–602, 2011.
- [46] L. Dagum and R. Menon, “OpenMP: an industry standard api for shared-memory programming,” Computational Science & Engineering, IEEE, vol. 5, no. 1, pp. 46–55, 1998.
- [47] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of mathematical functions. New York: Cambridge University Press, 2010.
- [48] A. Jeffrey and H.-H. Dai, Handbook of Mathematical Formulas and Integrals, 4th ed. Academic Press, 2008.

Doruk Tayli(S’13) received his B.Sc. degree in Electronics Engineering from Istanbul Technical University and his M.Sc. in degree in Communications Systems from Lund University, in 2010 and 2013, respectively. He is currently a Ph.D. student at Electromagnetic Theory Group, Department of Electrical and Information Technology at Lund University. His research interests are Physical Bounds, Small Antennas and Computational Electromagnetics. |

Miloslav Capek(SM’17) received his Ph.D. degree from the Czech Technical University in Prague, Czech Republic, in 2014. In 2017 he was appointed Associate Professor at the Department of Electromagnetic Field at the same university. He leads the development of the AToM (Antenna Toolbox for Matlab) package. His research interests are in the area of electromagnetic theory, electrically small antennas, numerical techniques, fractal geometry and optimization. He authored or co-authored over 70 journal and conference papers. Dr. Capek is member of Radioengineering Society, regional delegate of EurAAP, and Associate Editor of Radioengineering. |

Vit Losenicky received the M.Sc. degree in electrical engineering from the Czech Technical University in Prague, Czech Republic, in 2016. He is now working towards his Ph.D. degree in the area of electrically small antennas. |

Akrou Lamyae received the Dipl.-Ing. degree in networks and telecommunications from National School of Applied Sciences of Tetouan in 2012. Since 2014 she is working towards her Ph.D. degree in Electrical and Computer Engineering at the University of Coimbra. |

Lukas Jelinek received his Ph.D. degree from the Czech Technical University in Prague, Czech Republic, in 2006. In 2015 he was appointed Associate Professor at the Department of Electromagnetic Field at the same university. His research interests include wave propagation in complex media, general field theory, numerical techniques and optimization. |

Mats Gustafsson (SM’17) received the M.Sc. degree in Engineering Physics 1994, the Ph.D. degree in Electromagnetic Theory 2000, was appointed Docent 2005, and Professor of Electromagnetic Theory 2011, all from Lund University, Sweden. He co-founded the company Phase holographic imaging AB in 2004. His research interests are in scattering and antenna theory and inverse scattering and imaging. He has written over 90 peer reviewed journal papers and over 100 conference papers. Prof. Gustafsson received the IEEE Schelkunoff Transactions Prize Paper Award 2010 and Best Paper Awards at EuCAP 2007 and 2013. He served as an IEEE AP-S Distinguished Lecturer for 2013-15. |