Accurate and Efficient Evaluation of Characteristic Modes

Accurate and Efficient Evaluation of Characteristic Modes

Doruk Tayli,  Miloslav Capek,  Lamyae Akrou, Vit Losenicky, Lukas Jelinek, and Mats Gustafsson,  Manuscript received July 26, 2019; revised July 26, 2019. This work was supported by the Swedish Foundation for Strategic Research (SSF) under the program Applied Mathematics and the project Complex analysis and convex optimization for EM design, and by the Czech Science Foundation under project No. 15-10280Y.D. Tayli and M. Gustafsson are with the Department of Electrical and Information Technology, Lund University, 221 00 Lund, Sweden (e-mail: {doruk.tayli,mats.gustafsson} Capek, V. Losenicky and L. Jelinek are with the Department of Electromagnetic Field, Faculty of Electrical Engineering, Czech Technical University in Prague, Technicka 2, 166 27 Prague, Czech Republic (e-mail: {miloslav.capek,losenvit,lukas.jelinek} Akrou is with the Department of Electrical and Computer Engineering, Faculty of Sciences and Technology, University of Coimbra, Polo II, Pinhal de Marrocos, 3030-290 Coimbra, Portugal (e-mail:

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.

Antenna theory, numerical analysis, eigenvalues and eigenfunctions, electromagnetic theory, convergence of numerical methods.

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 dynamics111The 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


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


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]


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


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


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


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]


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


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


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


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


and the size of the matrix is , where


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


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


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


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 arithmetic222As 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]


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 times333Computations are done on a workstation with i7-3770 CPU @  GHz and  GB RAM, operating under Windows 7. of different matrices444Computation 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), .

Fig. 1: Convergence of the matrix  to the matrix  on the rectangular plate (Example R2) for different order of spherical modes and multiple electric sizes , with the highest spherical mode order . The superscript denotes the Frobenius norm. The convergence is computed with quadruple precision using the mpmath Python library [27].
Example Time to assemble matrices in IDA (s)
(see Table II)
TABLE I: Time to assemble matrices in IDA. Simulation setup for the examples in Table II, and (), matrix multiplication is performed with dgemm from the Intel MKL library [28].

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
Spherical shell S3
Fig. 9 S4
Rectangular plate R1
Fig. 10 R2
Helicopter H1
TABLE II: Summary of examples used throughout the paper, is the electrical size, is the number of basis functions (4), and is number of spherical modes calculated as (12). The order of the symmetric quadrature rule used to compute the non-singular integrals in (5) is  [34].
Number of properly calculated modes
(see Table II) (17) (19) (20) (24)
TABLE III: Comparison of the number of modes correctly found by the classical and the novel methods for examples listed in Table II. Columns summarize the radiation modes and columns summarize the CMs. Values in parentheses depicts the number of non-degenerated TM and TE modes found on spherical shell. The main outcome of the table, comparison of the CMs is highlighted by bold type.

Iii-a Radiation Modes

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


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


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


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 higher555Quantitatively, 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


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 


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


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


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


with expansion coefficients of CMs defined as


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


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 time666For 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)


and multiplying from the left with


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


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.

Fig. 2: The absolute values of the CMs of spherical shell with electrical size (Example S2). Data calculated with classical procedure (20) are compared with techniques from this paper, (21), (24), and with the analytical results valid for the spherical shell [12].
Fig. 3: The absolute values of the CMs of rectangular plate (Example R2). Since unknown analytical results, the multiprecision package Advanpix has been used instead to calculate the first  modes from impedance matrix in quadruple precision.
Example Time to calculate CMs (s)
(see Table II) (20) (24) (29)
TABLE IV: Comparison of computation time required by various methods capable to calculate first CMs. The calculations were done on Windows Server 2012 with 2Xeon E5-2665 CPU @  GHz and  GB RAM.
Fig. 4: Current density of the first characteristic mode of a helicopter at (Example H2), mesh grid has been taken from [32].

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.

Fig. 5: Similarity of numerically evaluated characteristic currents for a spherical shell of two different discretizations (Example S2 and S5) and the analytically known currents [12]. The coefficients  were calculated according to [12], top panel depicts results for the conventional procedure (20), bottom panel for the procedure from this paper (24).

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


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.

Fig. 6: Similarity of numerically and analytically evaluated characteristic far fields for a spherical shell of two different discretizations (Example S2 and S5) and analytically known far fields [12]. The coefficients  were calculated by (30) , see [12] for more details. Top panel depicts results for the conventional procedure (20), bottom panel for the procedure from this paper (24).

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.


[loop, controls, autoplay, buttonsize=0.5em, poster=28]1./TikZ_animate/OLD_merge030


[loop, controls, autoplay, buttonsize=0.5em, poster=28]1./TikZ_animate/NEW_merge030

Fig. 7: Comparison of the higher-order CMs of the rectangular plate (Example R2) with the most similar characteristic number, left panel: conventional procedure (20), right panel: procedure from this paper (24).

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.

Fig. 8: Left pane: the absolute values of the CMs of a spherical shell (Example S2) if only odd rows of the matrix  were kept. Right pane: the absolute values of the CMs of a rectangular plate (Example R2) if only even rows of the matrix  were kept.

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


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.


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]


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]


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

The spherical vector waves are [20, 13]


where are the radial function of order defined as

, (35a)
, (35b)
, (35c)

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


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


Appendix C Associated Legendre Polynomials

The associated Legendre functions are defined [48] as




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


The normalized associated Legendre function , is defined as follows


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


where , .

Appendix D Spherical Shell and Rectangular Plate

Meshes for the spherical shell of radius with and d-o-f are depicted in Fig. 9. The meshes for the rectangular plate of aspect ratio with , , and d-o-f are presented in Fig. 10.

Fig. 9: Spherical shell mesh with  triangles (left) and  triangles (right) with  (left) and  (right) RWG basis functions, respectively. The same mesh grids are used in [12] to make the results comparable.
Fig. 10: Rectangular plate mesh with , , and  triangles (from left to right) with , , and  RWG basis functions, respectively.

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., .

Fig. 11: Normalized eigenvalues of the matrix  of a spherical shell with electrical size discretized into  triangles (Example S2). Multiprecision package Advanpix [33] has been used for evaluation in quadruple precision. The number of well-determined modes is delimited by horizontal dashed lines.
Fig. 12: Normalized eigenvalues of the matrix  of rectangular plate (Example R2). Since the matrix  has only  rows, the number of modes is limited. The number of well-determined modes is delimited by horizontal dashed lines.


  • [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:
  • [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:
  • [28] Intel. (2017) Intel Math Kernel Library 2017 update 3. [Online]. Available:
  • [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:
  • [31] D. Tayli. (2017) IDA (Integrated Development toolset for Antennas). Lund University.
  • [32] Altair. (2016) FEKO. Altair. [Online]. Available:
  • [33] Advanpix. (2016) Multiprecision Computing Toolbox for MATLAB. [Online]. Available:
  • [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:
  • [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:
  • [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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description