-0.5

Y

0

0.5-1

-0.5
X

0

0.5

1.5

1

0.5

Z

θ(deg)
0 50 100 150 200

R
C
S
[d
B
sm

]

-25

-20

-15

-10

-5

0

5

p = 2

p = 4

p = 6

Mie

Method of Moments for Maxwell’s equa-
tions based on higher-order interpolatory
representation of geometry and currents

MAGNUS CARLSSON

Department of Signals and Systems
Chalmers University of Technology
Gothenburg, Sweden 2016





Master’s thesis EX066/2016

Method of Moments for Maxwell’s equations
based on higher-order interpolatory

representation of geometry and currents

MAGNUS CARLSSON

Department of Signals and Systems
Division of Signal processing and Biomedical engineering

Chalmers University of Technology
Gothenburg, Sweden 2016



Method of Moments for Maxwell’s equations based on higher-order interpolatory
representation of geometry and currents
MAGNUS CARLSSON

Copyright c© MAGNUS CARLSSON, 2016.

Supervisor and Examiner: Thomas Rylander, Department of Signals and Systems
Co-supervisor: Johan Winges

Master’s Thesis EX066/2016
Department of Signals and Systems
Division of Signal processing and Biomedical engineering
Chalmers University of Technology
412 96 Gothenburg
Telephone +46 31 772 1000

Cover: Upper left: Electric near field around two PEC sphere with incident hori-
zontal light.
Upper right: Electric near field around the PEC sphere with incident vertical light.
Down left: A curved surface element with a first-order Taylor approximation of the
surface.
Down right: A Radar Cross-Section spectrum as a function on incident polar angle
for a PEC sphere.

Typeset in LATEX
Gothenburg, Sweden 2016



Abstract

Electromagnetic scattering is a widely explored field in science and engineering,
where a vast number of applications can be found in the open litterature. The
Method of Moments (MoM) is an efficient numerical techniuqe for solving scattering
problems. This thesis explores the combination of a collocation scheme and the usage
of higher-order divergence-conforming basis functions for the currents induced on the
scatterer. The problem of a perfect electric conducting (PEC) sphere is analyzed and
the numerical results are in good agreement with analytical solutions. Exponential
convergence in the polynomial order p of the current expansion is achieved for p ≤ 6.
In addition, two adjacent PEC spheres are analyzed for an incident plane wave
with (i) the electric field polarized along the straight line between the centers of the
spheres and (ii) the electric field polarized perpendicular to the straight line between
the centers of the spheres.

i



ii



Acknowledgement

This thesis is the concluding part of the Master’s Program in Applied Physics at
Chalmers University of Technology. The thesis work has been conducted at the
department of Signals and Systems at Chalmers University of Technology during
autumn 2015 and beginning of spring 2016.

I would like to express my gratitude to everyone that in some way has contributed
to this thesis, and especially to the following persons:

Assoc. Prof. Thomas Rylander, my supervisor and examiner, for making this thesis
possible, for his great dedication in this work and for all the interesting and inspiring
discussions.

Assoc. Prof. Matthys Botha, for his contribution to the treatment of Singular
integrals and his great knowledge and guidance in the Method of Moments.

Johan Winges, for guidance, support, interesting discussions and for his contri-
bution to the programming part of this thesis.

Göteborg, April 2016
Magnus Carlsson

iii



iv



Notations

E Electric field (V/m)

D Electric flux density (C/m2)

H Magnetic field (A/m)

B Magnetic flux density (T)

J Electric current density (A/m2)

M Magnetic current density (V/m2)

A Magnetic vector potenial (Vs/m)

F Electric vector potenial (VF/m)

φe Electric potential (V)

φm Magnetic potential (A)

ρ Electric charge density (C/m3)

% Magnetic charge density (T/m)

ω Angular frequency (rad/s)

ε Permittivity (F/m)

µ Permeability (H/m)

j Imaginary unit

v



vi



Contents

1 Introduction 1
1.1 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Boundary conditions at interfaces between different media . . 2
1.2 Integral representations of the electromagnetic fields . . . . . . . . . . 3

1.2.1 Alternative representations (suitable for the near-field region) 6
1.2.2 Far-field representations . . . . . . . . . . . . . . . . . . . . . 6
1.2.3 Green’s functions . . . . . . . . . . . . . . . . . . . . . . . . . 7

2 Method of Moments 11
2.1 Integral equations for PEC scatterers . . . . . . . . . . . . . . . . . . 12
2.2 Representation of the geometry . . . . . . . . . . . . . . . . . . . . . 13

2.2.1 Discretization of a sphere . . . . . . . . . . . . . . . . . . . . 14
2.2.2 Reference element and higher-order mapping . . . . . . . . . . 15
2.2.3 Co- and contravariant basis vectors . . . . . . . . . . . . . . . 17
2.2.4 Surface integration . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3 Representation of the current density . . . . . . . . . . . . . . . . . . 18
2.3.1 Divergence conforming basis functions . . . . . . . . . . . . . 18

2.4 Weighting functions – Petrov-Galerkin’s method . . . . . . . . . . . . 20
2.5 Singular integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.5.1 Subdivision of the quadrilateral element . . . . . . . . . . . . 22
2.5.2 Mapping to suitable domain for product integration rules . . . 22
2.5.3 Integrands with weak singularity . . . . . . . . . . . . . . . . 25
2.5.4 Integrands with strong singularity . . . . . . . . . . . . . . . . 26

3 Results 27
3.1 Integration of known surface current . . . . . . . . . . . . . . . . . . 27

3.1.1 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Scattering from PEC sphere . . . . . . . . . . . . . . . . . . . . . . . 30
3.3 Scattering from two PEC spheres . . . . . . . . . . . . . . . . . . . . 33

4 Conclusions and future work 35
4.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Bibliography 36

A Far-field derivation for L and K 39

vii



B Projection point on tangent plane 43

C Strong true singularity treatment 45

D Analytic solutions for a test problem 49

viii



ix



Chapter 1

Introduction

Electromagnetic scattering caused by physical objects is a frequently occurring topic
in engineering and science. In radar [1,2] applications, the object may be an airplane
or a ship and significant scattering is often caused by metallic conductors that are
part of the construction. A possible interest could be to improve the stealth prop-
erties [3,4] of these objects. This could be achieved by investigating the radar cross
section (RCS) of the object. In other cases, such as light scattering from metallic
nanoparticles, plasmonic resonances [5, 6] occurs at certain frequencies. At these
frequencies, the scattering objects may change absorbing and transmitting spectra
dramatically.

To gain a deeper understanding of the previously mentioned topics, the theory
of electromagnetic scattering needs to be exploited. In practice, numerical methods
are needed since analytic solutions are very complicated or impossible to derive. The
Method of Moments (MoM) is an attractive technique used to solve electromagnetic
scattering problems in the frequency domain and it is the topic of this work. In the
following chapters, the theory of MoM is outlined and applied on a canonical test
case, namely the scattering of a perfect conducting sphere (PEC) sphere.

1.1 Maxwell’s equations

Maxwell’s equations describe the phenomena of classical electromagnetics, which in
the frequency domain are

∇×E = −M − jωB, (1.1)

∇×H = J + jωD, (1.2)

∇ ·D = ρ, (1.3)

∇ ·B = %, (1.4)

where we assume harmonic time dependence of the fields: E(r, t) = Re[E(r)ejωt],
H(r, t) = Re[H(r)ejωt]. The time dependence ejωt is supressed throughout this
thesis.

1



In equations (1.1) - (1.4), E is the electric field, B is the magnetic flux density,
H is the magnetic field, D is the displacement field, J is the electric current density
and ρ is the electric charge density. The magnetic current density M and charge
density % are not physically realizable quantities but are useful mathematical tools
in solving radiation and scattering problems.

The continuity equation for electric and magnetic charge are,

jωρ+∇ · J = 0, (1.5)

jω%+∇ ·M = 0. (1.6)

This follows from taking the divergence of equations (1.1)-(1.2) and combining with
equations (1.3)-(1.4).

In this thesis, we restrict our attention to linear and isotropic materials. For
such materials, the consitutive relations are

B = µH , (1.7)

D = εE, (1.8)

where ε is the electric permittivity and µ is the magnetic permeability for the specific
material that is considered.

1.1.1 Boundary conditions at interfaces between different
media

Consider two regions, V1 and V2, separated by the surface interface S. Let the
normal vector n̂2 be perpendicular to S such that it points into V1. On the interface
S, we have the boundary conditions

n̂2 × (E1 −E2) = −M s, (1.9)

n̂2 × (H1 −H2) = J s, (1.10)

n̂2 · (D1 −D2) = ρs, (1.11)

n̂2 · (B1 −B2) = %s, (1.12)

where the surface current and charge densities coincide with the interface S between
V1 and V2. Equations (1.9) and (1.10) state the discontinuity in the tangential parts
of H and E at the interface between the two regions are equal to the surface
current densities at the interface. Equations (1.11) and (1.12) state similary that
the discontinuity between the normal components of D and B are equal to the
surface charge densities at the interface.

If V2 is a perfect electric conductor (PEC) and V1 is a dielectric, the boundary
conditions reduce to

2



n̂2 ×E1 = 0, (1.13)

n̂2 ×H1 = J s, (1.14)

n̂2 ·D1 = ρs, (1.15)

n̂2 ·B1 = 0, (1.16)

since all the fields inside the PEC are zero at non-zero frequency. The PEC boundary
conditions are a very good approximation for metallic surfaces in the microwave
regime.

1.2 Integral representations of the electromagnetic

fields

In the following section, we rewrite Maxwell’s equaitons (1.1) - (1.4) as integral
equations, where the scattered fields Es and Hs may be expressed as functions of
the induced surface currents J s andM s on the scatterer. This is done by introducing
so-called Green’s functions.

We illustrate this procedure in the context of electrostatics [7]. Consider the
inhomogeneus differential scalar Poisson’s equation

∇2φ = − ρ
ε0
, (1.17)

for the electric potential φ. For a system of n point charges qi located at r′i, with
i = 1, 2, ..., n, the charge density can be expressed as

ρ =
n∑
i=1

qiδ
3(r − r′i). (1.18)

The solution of Poisson’s equation (1.17) at r in free space can be constructed by
superposition:

φ(r) =
n∑
i=1

qi
4πε0|r − r′i|

. (1.19)

For a continuum charge density ρ(r′), the sum is replaced by an intergral

φ(r) =

∫
V

ρ(r′)

4πε0|r − r′|
dV ′. (1.20)

It is possible to show that (1.20) and (1.17) are equivalent by applying the Laplace
operator on both sides in (1.20). The resulting integrand is singular, but by noting
that [7]

∇2

(
1

|r − r′|

)
= −4πδ3(r − r′), (1.21)

3



we find the expressions are indeed equivalent. For this specific problem, we identify
form (1.20) that the Green’s function is

G(r, r′) =
1

4ε0π|r − r′|
, (1.22)

which represents the potential or field at r produced by a point charge at r′. Hence,
the differential equation formulation (1.17) could be reformulated as an integral
equation

φ(r) =

∫
V

G(r, r′)ρ(r′) dV ′, (1.23)

where the objective is to compute the source distribution ρ given the field φ. This
procedure can be generalized to include other physical situations. Consider a differ-
ential equation

D [f ] = s, (1.24)

where D is a differential operator, f is a scalar field (or a component of a vector field)
and s is a scalar source distribution (or a component of a vector source distribution).
The Green’s function G(r, r′) represents the value of f at r due to a point source
at r′ and it satisfies

D [G(r, r′)] = δ3(r − r′). (1.25)

Multiplying (1.25) with s(r′) and integrating with respect to r′ yields∫
D [G(r, r′)s(r′)] dV ′ =

∫
δ3(r − r′)s(r′) dV ′. (1.26)

The right-hand side reduces to s(r) due the nature of the Dirac distribution δ3(r−r′)
and the left-hand side can be rewritten as

D

[∫
G(r, r′)s(r′) dV ′

]
, (1.27)

due to the linearity of the differential operator. By comparison with (1.24), we can
identify the integral formulation

f(r) =

∫
G(r, r′)s(r′) dV ′, (1.28)

of the initial differential formulation. It is clear that the integration domain in (1.28)
reduces to the region where sources exist. This technique is thus very powerful for
localized source distributions in an otherwise source-free region. This is often the
case in scattering problems, where a relatively small metallic object scatters a field
into free space. For conductors, the source currents are generally restricted to the
surface of the conductor [8] and, thus, the sources are located only on the surface of
the scatterer. In the following, we assume that there is only currents on the surface
of the geometry and we denote the surface current densities by J and M to simplify
the notation. The fields are computed at the obervation point r and the sources are
located at r′.

4



For homogeneous space with permittivity ε and permeability µ, the above proce-
dure can be applied to find the integral equations for the scattered fields produced
by the surface currents J and M , where a detailed derivation can be found in [9].
The scattered fields Es and Hs may be expressed as

Es = −jωA−∇φe −
1

ε
∇× F , (1.29)

Hs = −jωF −∇φm +
1

µ
∇×A, (1.30)

where A and F are the magnetic and electric vector potentials, φe and φm are
the electric and magnetic scalar potential, respectively. The expressions for these
quantities are

A(r) = µ

∫
S

J(r′)G(r, r′) dS ′, (1.31)

F (r) = ε

∫
S

M (r′)G(r, r′) dS ′, (1.32)

φe(r) =
1

ε

∫
S

ρ(r′)G(r, r′) dS ′, (1.33)

φm(r) =
1

µ

∫
S

%(r′)G(r, r′) dS ′, (1.34)

where G(r, r′) is the dynamic Green’s function in 3D for a homogeneous region that
is unbounded and the integration is performed over the surface S of the conductor.
The Green’s function is

G(r, r′) =
e−jkr

4πr
, (1.35)

where r = |r − r′| and k = ω
√
εµ. This Green’s function is derived from the

inhomogeneous Helmholtz equation

−
(
∇2f(r) + k2f(r)

)
= s(r). (1.36)

The scattered fields may be expressed in terms of the linear operators L and K
as sole functions of J and M , by inserting (1.31) - (1.34) in (1.29) - (1.30) and using
the expressions for the continuity of charge and current density (1.5) - (1.6), which
yields

Es(r) = −jωµ(LJ)(r)− (KM )(r), (1.37)

Hs(r) = −jωε(LM)(r) + (KJ)(r), (1.38)

where the linear operators L and K are defined as

(LX)(r) =

[
1 +

1

k2
∇∇ ·

] ∫
S

G(r, r′)X(r′) dS ′, (1.39)

(KX)(r) = ∇×
∫
S

G(r, r′)X(r′) dS ′. (1.40)

The integral equations (1.37) and (1.38) are general expressions for the scattered
fields. However, the scattered fields are usually divided into the near-field and far-
field regions, where terms of different order O(rn) are dominating.

5



1.2.1 Alternative representations (suitable for the near-field
region)

The near-field region comprises the fields produced close to the sources. For field
points that are close to the source distribution, alternative expressions of (1.39) and
(1.40) may be used. The derivative operators are moved inside the integral signs,
which yield

(LX)(r) =

∫
S

G(r, r′)X(r′) dS ′ +
1

k2

∫
S

∇∇ · [G(r, r′)X(r′)] dS ′, (1.41)

(KX)(r) =

∫
S

∇× [G(r, r′)X(r′)] dS ′, (1.42)

where all derivatives are with respect to r. Since X(r′) only depends on r′, the
derivatives only affects the Green’s function. This yields terms that scale as 1/r, 1/r2

and 1/r3. Thus, the Green’s function exhibits a singularity as the distance between
source point and field tends to zero. Special care is neccesary for numerical integra-
tion of these terms as the results otherwise may become very inaccurate. Usually,
the singularities are denoted as weak singularities (1/r), strong singularities (1/r2)
and hyper singularities (1/r3). Special care is also taken if r and r′ coincide exactly
(true singularity) or are very close (near singularity).

Integration schemes designed to handle the near-field expressions can be found
in the open literature for weak singularities [10] and for strong singularities [11].
Hyper singularities can be reduced to strong singularities by transferring one of the
derivatives to the source. The second integral in (1.43) possess a hyper singularity,
which after integration by parts [9] yields the operators

(LX)(r) =

∫
S

G(r, r′)X(r′) dS ′ +
1

k2

∫
S

∇G(r, r′)∇′ ·X(r′) dS ′, (1.43)

(KX)(r) =

∫
S

[∇G(r, r′)]×X(r′) dS ′, (1.44)

and the integration schemes for weak and strong singularities suffice. Here, the
differation formula ∇× (GX) = ∇G×X ′+G∇×X ′ = ∇G×X ′ is used for (1.44).

1.2.2 Far-field representations

When r is located sufficiently far away from r′, we approximate r = |r − r′| as [7]

r =

{
r − r̂ · r′ for phase variations

r for amplitude variations
, (1.45)

given the far-field conditions r � d, r � λ and r � kd2, where d = max|r′| is the
maximum extension of the scatterer and λ = 2π/k. Consider the expressions for the
scattered fields in (1.37) and (1.38). When acting upon the Green’s function with
the differential operators, it produces terms in the fields that vary according to 1/r,
1/r2 and 1/r3. In the far-field, only the fields that vary as 1/r have a significant

6



amplitude. In addition, the electromagnetic wave is transverse i.e. there are no field
components along the direction of propagation. The far-field expressions for the L-
and K-operators then reduce to

(LX)(r) = −r̂ × r̂ × P (r̂)
e−jkr

kr
, (1.46)

(KX)(r) = −jkr̂ × P (r̂)
e−jkr

kr
, (1.47)

where

P (r̂) =
k

4π

∫
S

e−jkr̂ ·k̂X(r′) dS ′, (1.48)

is the spatial Fourier transform of the current density X evaluated in the coordinate
kr̂. A detailed derivation of the far-field expressions is given in Appendix A.

1.2.3 Green’s functions

The physical interpretation of the Green’s function is the field at r caused by a unit
point source at at r′ [12]. Below, we provide a derivation of the dynamic Green’s
function and introduce the so-called dyadic Green’s function.

More generally, it entered as an integral solution (1.28) of an initial differential
equation (1.24). The derivation of the Green’s function is based on that the solution
can be expresses as the superposition of the contribution from many point sources,
all weighted by G(r, r′).

Scalar Green’s function

The scalar Green’s function can be derived from

(∇2 + k2)G(r, r′) = −δ3(r − r′). (1.49)

where G(r, r′) is a solution to the inhomogeneous Helmholtz equation (1.36). Since
we have assumed homogeneous space, G(r, r′) may only depend on the distance
between the source and observation points r = |r − r′|. The Laplacian in spherical
coordinates for a scalar field f is given by

∇2f =
1

r2

∂

∂r

(
r2∂f

∂r

)
+

1

r2 sin θ

∂

∂θ

(
sin θ

∂f

∂θ

)
+

1

r2 sin2 θ

∂2f

∂ϕ

2

. (1.50)

It is then possible to rewrite (1.49) as

1

r2

d

dr

(
r2 dG(r)

dr

)
+ k2G(r) = 0, r > 0. (1.51)

The solution of (1.51) may be found in [13] for example and has the form

G(r) = A
e−jkr

r
, (1.52)

7



where A is an arbitrary constant and we have assumed only outward traveling waves
from the source. The constant may be determined by substituting (1.52) in (1.49)
and integrate over a sphere with radius r0, where we let r0 tend to zero. The different
terms evaluate to∫

r<r0

∇2G dV =

∮
r=r0

∇G · n̂dS = 4πr2
0r̂ · ∇G

= −4πAr2
0

(
e−jkr0

r2
0

+ jk
e−jkr0

r0

)
∫
r<r0

k2G dV = 4πk2A

[
− 1

jk
r0e
−jkr0 +

1

k2
(e−jkr0 − 1)

]
∫
r<r0

−δ3(r) dV = −1.

By taking the limit r0 → 0, we find that A = 1/4π and we finally get the scalar
Green’s function

G(r, r′) =
e−jkr

4πr
. (1.53)

Dyadic Green’s functions

The expressions for the scattered fields (1.43) and (1.44) may be used when the
distance between the source and field point is sufficiently large. Since the derivatives
are only acting on the unprimed coordinates they only act on the scalar Green’s
function G(r, r′). The result may be expressed in terms of dyadic Green’s functions
according to

¯̄GL(r, r′) =

[
¯̄I +

1

k2
∇∇

]
G(r, r′), (1.54)

¯̄GK(r, r′) = ∇×G(r, r′). (1.55)

In this thesis, we perform the integrations in Cartesian coordinates. The dyadic
Green’s functions can then be expressed in matrix form and are given by

¯̄GL(r, r′) =


k2 + ∂2

∂x2
∂2

∂x∂y
∂2

∂x∂z

∂2

∂y∂x
k2 + ∂2

∂y2
∂2

∂y∂z

∂2

∂z∂x
∂2

∂z∂y
k2 + ∂2

∂x2

 G(r, r′)

k2
, (1.56)

¯̄GK(r, r′) =


0 − ∂

∂z
∂
∂y

∂
∂z

0 − ∂
∂x

− ∂
∂y

∂
∂x

0

G(r, r′). (1.57)

8



The Green’s function is then moved inside the matrix and the differentiations are
calculated analytically. This operation is straightforward but algebraically tedious.
The results are given below

¯̄GL(r, r′) =


k2 + Fxx Fxy Fxz

Fyx k2 + Fyy Fyz

Fzx Fzy k2 + Fzz

 G(r, r′)

k2
, (1.58)

¯̄GK(r, r′) =


0 −Fz Fy

Fz 0 −Fx
−Fy Fx 0

G(r, r′), (1.59)

where we use the following functions

Fm =
(ξm − ξ′m)(−1− jkR)

R2
, (1.60)

Fmn =
(ξm − ξ′m)(ξn − ξ′n)(3(1 + jkR)− k2R2)

R4
, (1.61)

Fmm =
3(ξm − ξ′m)2

R4
+

3jk(ξm − ξ′m)2

R3
− 1 + k2(ξm − ξ′m)2

R2
− jk

R
, (1.62)

where ξm, ξn = x, y or z. Inserting the expressions for the dyadic Green’s function
in (1.37) and (1.38) yield a compact notation

(LX)(r) =

∫
S

¯̄GL(r, r′)X(r′) dS ′, (1.63)

(KX)(r) =

∫
S

¯̄GK(r, r′)X(r′) dS ′. (1.64)

9



10



Chapter 2

Method of Moments

The Method of Moments (MoM) is a technique to convert the integral equations
decribed in Chapter 1 into a system of linear equations. It bears strong resemblance
to the Finite Element Method (FEM), which is used to solve differential equations.
To construct the set of linear equations, we use the weighted residual method, a
general recipe for both the FEM and MoM.

For the MoM, the vector integral equations (1.37) and (1.38) have the form

L[X] = V , (2.1)

where L denotes a linear integral operator, V is a given excitation function, and X
is the unknown current distribution to be solved on the surface S. The surface S is
discretized into elements, which for example may be triangles or quadrilaterals. The
solution X is expanded in a linearly independent vector basis functions according
to

X(r) ≈
N∑
j=1

XjBj(r), (2.2)

where Xj are the unknown coefficients to be determined. The basis functions Bj(r)
are low-order polynomials that are non-zero only in a few adjecent elements. We
form the residual

r = L[X]− V , (2.3)

which we want to minimize. This may be done by introducing weighting functions
wi (as many as there are unknown coefficients) for weighting the residual (2.3). By
enforcing the average of the weighted residual to be zero, i.e.∫

S

wi · r dS = 0, (2.4)

and inserting the expansion (2.2), we can solve for the unknown coefficients Xj

N∑
j=1

Xj

∫
S

wi · L[Bj(r)] dS =

∫
S

wi · V dS. (2.5)

Of particular importance is the choice of weighting functions wi. A common ap-
proach is to choose the weighting functions to be the same as the basis functions

11



wi = Bi, which is referred to as Galerkin’s method. Another approach is the so-
called point matching method. It corresponds to taking the weighting functions as
delta functions wi = δ3(r− ri)ξ, where ξ is a tangential vector to the surface. This
choice enforces the residual to be zero at the matching points. The point matching
method is exploited in this thesis.

2.1 Integral equations for PEC scatterers

The scattered fields are given by (1.37) and (1.38). We impose the boundary con-
ditions on the PEC according to (1.13) and (1.14), where E and H are the total
tangential fields on the surface. The total field is the superposition of the incoming
and scattered field, where the scattered field is caused by the induced surface current
on the scatterer

n̂ × (Ei +Es(J)) = 0, (2.6)

n̂ × (H i +Hs(J)) = J . (2.7)

Equation (2.6) states that the total tangential electric field on the conductor must
vanish on the surface. Likewise, (2.7) states that the total tangential magnetic field
is related to the electric current. These boundary conditions may be enforced by
scalar multiplying the total field with a tangential weighting function and integrate
over the surface S according to∫

S

wi · (Ei +Es(J)) dS = 0, (2.8)∫
S

wi · (H i +Hs(J)) dS =

∫
S

wi · J dS. (2.9)

We may then insert the expressions for the scattered fields in terms of the operators
L and K ∫

S

wi · (−jωµL(J)) dS = −
∫
S

wi ·Ei dS, (2.10)∫
S

wi · (K(J)) dS −
∫
S

wi · J dS = −
∫
S

wi ·H i dS, (2.11)

where we have M = 0, since we consider a PEC scatterer. Equation (2.10) is known
as the electric field integral equation (EFIE) and (2.11) is known as the magnetic
field integral equation (MFIE). Once J is determined by (2.10) or (2.11), the fields
and other observable quantities associated with electromagnetic scattering can be
computed.

We use the expansion in (2.2) to discretize the current according to the MoM.
Inserting this expression in (2.10) and (2.11) yields a system of linear equations

Zi = b, (2.12)

12



where the elements of the impedance matrices are given as

ZEFIE
ij =

∫
S

wi · (−jωµL(Bj)) dS, (2.13)

ZMFIE
ij =

∫
S

wi · (K(Bj)) dS −
∫
S

wi ·Bj dS, (2.14)

and the rows of the excitation vectors as

bEFIEi = −
∫
S

wi ·Ei dS, (2.15)

bMFIE
i = −

∫
S

wi ·H i dS. (2.16)

Once Z and b are assembled, the vector i of unknown coefficients may be obtained
with matrix inversion.

On closed conductors, it is a well known problem that the EFIE and MFIE cannot
produce unique solutions for all frequencies [8,9]. This corresponds to homogeneous
solutions to these equations that satisfy the boundary conditions in the absense of
an incoming field. These solutions are interior modes to the conductor and they
produce a singular matrix Z at their resonance frequencies. A technique to avoid
this problem is through a linear combination of the EFIE and the MFIE equations

ZCFIE = αZEFIE + (1− α)Z0Z
MFIE, (2.17)

bCFIE = αbEFIE + (1− α)Z0b
MFIE (2.18)

where α ∈ [0, 1] is a weighting parameter and this new scheme is called combined
field integral equation (CFIE). Since the EFIE and MFIE have different internal
resonance modes, α may be chosen to eliminate the singular behavior due to the
internal resonances [8, 9].

The incoming fields Ei and H i are assumed to originate from sources very far
away. Thus, the incident field can be approximated as a plane wave. The electric
field of a plane wave propagating in k̂ direction has the form

Ei = E0e
−jk·r, k = k̂k, (2.19)

where r = x̂x + ŷy + ẑz is the coordinate vector and E0 is a complex amplitude
with a direction that lies in the plane perpendiculary to the direction of propagation
k̂. Under the plane wave assumption, the magnetic field associated with the electric
field (2.19) is obtained by

H i =
1

Z
k̂ ×Ei, (2.20)

where Z =
√
µ/ε is the impedance of the surrounding medium.

2.2 Representation of the geometry

In this thesis, we assume the geometry of interest has no sharp edges or corners. It
is possible to include sharp edges or corners when formulating the EFIE and the
MFIE equations but it is outside the scope of this thesis.

13



The surface of the geometry is discretized into curved quadrilateral cells, where
the cells are also referred to as patches in the following. A quadratic approximation
with 3× 3 interpolatory nodes of a curved surface element Se in shown in Fig. 2.1.

-

6

�
�
�	

z

x

y

r r r
r r r

r r r
Se

Figure 2.1: A curvelinear cell Se in x-y-z space.

Each quadrilateral cell is described by a set of interpolatory nodes that coincide
with the original geometry’s surface. Another possibility to resolve the geometry is to
use triangles instead of quadrilaterals. However, quadrilaterals are more suitable for
numerical integration by means of product rules based on one-dimensional Gaussian
quadrature rules.

2.2.1 Discretization of a sphere

A spherical scatterer of a homogeneous, isotropic and linear medium allows for an-
alytical expressions on closed form that yield the scattered field given an incident
plane wave. Thus, it is a simple but very interesting test case for scattering com-
putations in electromagnetics. There are various ways to discretize the surface of a
sphere and, in this thesis, we discretize the surface of a cube and map it onto the
surface the sphere. Each face of the cube is discretized into rectangles, which yield
curved patches on the sphere. Of course, the discretization of the cube influences
the distribution of curved patches on the sphere.

It is easy and straightforward to discretize each face of the cube by a uniform and
structured mesh of squares. Consider a cube centered at the origin. The coordinate
vector of a point on the surface of the cube can be denoted rc. The corresponding
coordinate rs on the surface of a sphere with radius a can be obtained with the
transformation

rs = a
rc

|rc|
. (2.21)

An example of this discretization is shown Figs. 2.2(a) and 2.2(b). Each face of the
cube is discretized into 49 uniform sized squares and are then mapped to the sphere.

It is clear from Fig. 2.2(b) that the elements on the sphere are not uniform in
size. This is because the patches near the center of each face are more streched than
the patches along the edges. This may be a problem in the context of convergence
studies with respect to the cell size. For such studies, it is desireable to preserve the
uniformity in cell size on the spherical surface as the mesh is refined.

14



(a) Discretization of cube. (b) Discretization of sphere.

Figure 2.2: Mapped discretization of cube to surface of sphere.

A method to alleviate this problem is to use an area preserving mapping from
the cube to the sphere [14], by the means of inverse Lambert azimuthal equal area
projection. An example mesh for the sphere, based on 49 cells on each face of the
cube, is shown in Fig. 2.3.

Figure 2.3: Area preserving mapping of cube to surface of sphere.

2.2.2 Reference element and higher-order mapping

An efficient framework to both describe and perform surface integration over an
arbitrary surface element Se in R3 Euclidian space is through the concept of a
reference element. Consider the coordinate vector r in physical space (x, y, z) as a
function of the reference element coordinates (u, v), which is expressed as

r(u, v) = x(u, v)x̂ + y(u, v)ŷ + z(u, v)ẑ . (2.22)

15



s

s

s

s

s

s

s

s

s

-

6

v

u
u = 1

v = 1

u = −1

v = −1

Figure 2.4: Integration limits of [−1, 1]2 cell.

To construct the mapping r(u, v), we use Lagrangian interpolation of the physical
node coordinates rij according to

r(u, v) =
M∑
i,j

rijSij(u, v). (2.23)

Here, we use a 2D Lagrangian polynomial Sij(u, v) by forming the product of two
1D polynomials PM

i (u) and PM
j (v) as

Sij(u, v) = PM
i (u)PM

j (v), i, j = 1, 2 . . . ,M, (2.24)

where the 1D Lagrangian polynomial of order M − 1 is given by

PM
i (u) =

M∏
j=1,j 6=i

u− uj
ui − uj

, (2.25)

where −1 ≤ u ≤ 1 with uj uniformly distributed on the interval [−1, 1] and

PM
i (uj) =

{
1 if i = j

0 if i 6= j
. (2.26)

Hence, we need to provide M × M nodes from each cell. It is also possible to
calculate the derivatives {∂x

∂u
, ∂x
∂v
, ∂y
∂u
, ∂y
∂v
, ∂z
∂u
, ∂z
∂v
} from (2.23), which are used in the

Jacobian matrix

J =

∂x∂u ∂y
∂u

∂z
∂u

∂x
∂v

∂y
∂v

∂z
∂v

 . (2.27)

16



2.2.3 Co- and contravariant basis vectors

The covariant basis vectors a and b are tangential to the surface patch and they
may be expressed as a function of u and v as

a = ∂r
∂u

= ∂x
∂u
x̂ + ∂y

∂u
ŷ + ∂z

∂u
ẑ ,

b = ∂r
∂v

= ∂x
∂v
x̂ + ∂y

∂v
ŷ + ∂z

∂v
ẑ ,

(2.28)

where the derivatives are calculated from (2.23). If v is held constant and u is varied
in (2.22), we get a curve on the surface and a is tangential to this curve. Similary,
b is tangential to a curve where u is constant. The covariant basis vectors a and b
are a suitable choice for representing the surface current in the MoM.

Next, we define the contravariant basis vectors by

α = ∇u = ∂u
∂x
x̂ + ∂u

∂y
ŷ + ∂u

∂z
ẑ ,

β = ∇v = ∂v
∂z
x̂ + ∂v

∂y
ŷ + ∂v

∂z
ẑ .

(2.29)

The contravariant basis vectors are also tangential to the curved surface and, in
addition, they are normal to constant-u or constant-v curves on the surface. It can
be shown that (2.28) and (2.29) fullfill the orthogonality relations

a ·α = 1, b ·α = 0,

a · β = 0, b · β = 1.
(2.30)

2.2.4 Surface integration

The surface integrals described in Chapter 1 are expressed in terms of an integration
domain that extends over the entire surface S. To make use of the discretization
in Section 2.2, an arbitrary integral of a function f(r) over the surface S may be
divided into a sum of integrals over the surface elements Se according to∫

S

f(r) dS =
∑
e

∫
Se
f(r) dS. (2.31)

To integrate a quantity over a physical surface patch Se, we use the mapping de-
scribed by (2.23). If we make a change of variables to the reference space, we have
to include the determinant of the Jacobian matrix of the mapping in the integration.
It can be shown [15] that |a× b| corresponds to det J for a 2D surface in 3D space
and we use the notation D(u, v) = |a× b| in the following. Thus, a surface integral
of a function f(r) may be transformed to the reference cell in u− v space according
to ∫

Se
f(r) dS =

∫ 1

−1

∫ 1

−1

f(r(u, v))D(u, v)dudv. (2.32)

17



The integral in the right-hand side of (2.32) can be calculated with numerical quadra-
ture. Consider an integral of a function f(x) over the interval [−1, 1]. A numerical
quadrature rule of order N approximates the integral according to∫ 1

−1

f(x) dx ≈
N∑
p=1

f(xp)wp, (2.33)

where wp are weights and xp are quadrature points associated with a specific quadra-
ture rule. This can easily be extended to a 2D integral according to∫ 1

−1

∫ 1

−1

f(u, v) dudv ≈
N∑

p=1,q=1

f(up, vq)wpwq. (2.34)

In this thesis, we consider the Gauss-Legendre [16] and Gauss-Lobatto [17]
quadrature schemes. It can be shown that Gauss-Legendre quadrature with N
points provides the exact result if f(x) is a polynomial of degree 2N − 1 or less.
Similary, Gauss-Lobatto quadrature integrates exactly for polynomials of degree
2N − 3 or less. If f(x) is not a polynomial function but infinitively differentiable at
a real or complex number a, it can be expanded in Taylor series according to

f(x) = f(a) +
f ′(a)

1!
(x− a) +

f ′′(a)

2!
(x− a)2 +

f ′′′(a)

3!
(x− a)3 + . . . (2.35)

where Gauss-Legendre quadrature evaluates the first 2N−1 terms exactly. Similary,
Gauss-Lobatto quadrature evaluates the first 2N − 3 terms exactly.

2.3 Representation of the current density

Next, we construct a representation of the current distribution flowing on the surface.
One important physical consideration is the continuity equation for the current and
charge distribution (1.5). This equation must also be fulfilled when moving over
a cell boundary to an adjacent cell. This may be achieved by representing the
current distribution in terms of divergence-conforming basis functions. A divergence-
conforming basis function enforces continuity for the normal component of the vector
field across cell egdes [15]. A divergence-conforming representation of the surface
current yields a finite charge distribution for the entire surface. However, the charge
distribution may be discontinuous along the edges of the cells.

A particular set of basis functions that satisfy these conditions are so-called
Rao-Wilton-Glisson (RWG) functions [18], which are the lowest-order divergence-
conforming basis functions. The RWG basis functions for triangle and quadrilateral
cells can be found in [15].

2.3.1 Divergence conforming basis functions

The currents J and M are expanded according to (2.2). We may express the
divergence-conforming basis function Bj in terms of basis functions formulated on
the reference shown in Fig. 2.4, and then map those to the physical element in the

18



3D space (x, y, z). On the reference element, we construct the basis functions in
terms of Lagrange polynomials as

B̃
U

ij = Lpi (u)Gp−1
j (v)û, (2.36)

B̃
V

ij = Gp−1
i (u)Lpj(v)v̂, (2.37)

where (i, j) are local indices, Lp(u) is a Lagrange polynomial of order p with in-
terpolatory points that coincide with the quadrature points of the Gauss-Lobatto
integration scheme. Similary, Gp−1(u) is a Lagrange polynomial of order p− 1 with
interpolatory points that coincide with the quadrature points of the Gauss-Legendre
integration scheme. An example of the interpolatory degrees of freedom is shown

in Figs. 2.5(a) for B̃
U

ij and 2.5(b) for B̃
V

ij , where the order of the basis functions is
p = 2.

r

r

r

r

r

r

-

-

-

-

-

-

-

6

v

u

(a) B̃
U

ij = L3
i (u)G2

j (v)û

r r

r r

r r

6 6

6 6

6 6

-

6

v

u

(b) B̃
V

ij = G2
i (u)L3

j (v)v̂

Figure 2.5: Interpolatory degrees of freedom for B̃
U

ij (a) and B̃
V

ij (b) shown on the reference
element. Here, quadratic basis functions (p = 2) are shown.

We then need to map these functions to a specific element e. To account for the
curvature of element e, it can be shown [8,15] that the relation between a divergence-

conforming basis function B̃ij on the reference element and a divergence-conforming
basis function Bij,e on the physical element is given by

Bij,e =
[Je]T

De
B̃ij, (2.38)

where Je is given by Eq. (2.27) evaluated for element e. Equations (2.36) and (2.37)
yield the basis functions expressed in physical space according to

BU
ij,e =

1

De(u, v)
Lpi (u)Gp−1

j (v)ae, (2.39)

BV
ij,e =

1

De(u, v)
Gp−1
i (u)Lpj(v)be, (2.40)

19



where ae and be are the covariant basis vectors evaluated for element e accodring
to the expressions in Section 2.2.3. The global current X may then be expressed as
a sum of basis functions multiplied with unknown coefficients

X(r) =
∑
e

(∑
i,j

XU
ij,eB

U
ij,e +

∑
i,j

XV
ij,eB

V
ij,e

)
. (2.41)

Each basis function associated with a corresponding degree of freedom that is interior
the the element, i.e. does not coincide with its edges, is non-zero only on the element
where it is defined. In case of an edge basis function, it is non-zero on the two
adjacent elements that share the edge associated with the basis function. Thus, it
is useful to define a global index ((i, j)local, e) 7→ jglobal to relate the indices (i, j)local

that are local to element e to the corresponding global basis function.
By expanding the current according to (2.41), one may simplify the calculation

of integrals containing the dyadic Green’s functions dramatically. Inserting (2.41)
in either (1.63) or (1.64) yields

∑
e

∫
Se

¯̄G(r, r′)

(∑
i,j

XU
ij,eB

U
ij,e(r

′) +
∑
i,j

XV
ij,eB

V
ij,e(r

′)

)
dS ′. (2.42)

The integral over a patch Se is mapped to the reference element and integrated
numerically as described Section 2.2.4. With Eqs. (2.39) and (2.40) we arrive at∑

e

∑
p,q

[
¯̄G(r, r′(up, vq))

(∑
i,j

XU
ij,eL

p
i (up)G

p−1
j (vq)ae(up, vq)+∑

i,j

XV
ij,eG

p−1
i (up)L

p
j(vq)be(up, vq)

)]
wpwq

(2.43)

If we choose a quadrature scheme with quadrature points that coincide with the
interpolation point for the vector basis function, i.e. a suitable product rule of Gauss-
Legendre and Gauss-Lobatto quadrature, Eq.(2.43) simplifies to∑

e

(∑
i,j

XU
ij,e

¯̄G(r, r′(ui, vj))ae(ui, vj)wiwj+∑
i,j

XV
ij,e

¯̄G(r, r′(ui, vj))be(ui, vj)wiwj

)
,

(2.44)

since basis functions satisfy the relations

Bij,e =

{
1 if i = p and j = q

0 otherwise
. (2.45)

2.4 Weighting functions – Petrov-Galerkin’s method

We use the Petrov-Galerkin’s method and define the weighting functions

wi = Cδ2(r − ri)ξ (2.46)

20



where i is a global degree of freedom index, C is a normalizing constant to be
determined and ξ is a contravariant basis vector to the surface, which is given by
(2.29). The coordinate vector ri is chosen to be collocated with the intepolation
point for the corresponding degree of freedom Bi. This way, the residual (2.3) is
enforced to be zero at discrete points collocated with the degrees of freedom.

To find the constant C, we integrate the scalar part of the weighting function,
i.e. wi = Cδ2(r − ri), over the surface. It is desirable to set this integral to unity∫

S

wi dS = 1, (2.47)

for normalization purposes. Inserting the scalar part of (2.46) in (2.47), rewriting
the integral to a sum over the elements e and transforming to the reference element
yields ∑

e

∫ 1

−1

∫ 1

−1

Cδ2(r(u, v)− r(ui, vi))D
e(u, v)dudv. (2.48)

This evaluates to

CDe(ui, vi) = 1⇒ C =
1

De(ui, vi)
. (2.49)

If the degree of freedom is located at an edge between two cells, we choose the average
of two single weighting functions defined on the adjecent cells e and e′ according to

wi =
1

2
δ2(r − ri)

(
1

De
ξe +

1

De′
ξe

′
)
. (2.50)

By choosing the weighting functions wi according to the Petrov-Galerkin’s scheme,
the outer integral in the Eqs. (2.10) and (2.11) is reduced such that its integrand is
evaluated at a single point.

2.5 Singular integrals

When the field point and the source point are sufficiently separated, Gaussian
quadrature can be used directly for integration. However, the assembling proce-
dure frequently involves field points that are located on the surface. This is known
as the true-singularity case, since the integrand in the EFIE or MFIE equations
becomes singular (weak or strong) due the Green’s function.

Another case is when field points are located very close to the surface. In this
case, the integrand becomes nearly singular (weak or strong) and this is referred to
as the near-singularity case.

One technique for numerically evaluating integrands for these cases is to trans-
form the integration domain so that the Jacobian of the transformation cancels the
singularity produced by the factor 1/rn of the Green’s function. There exists many
transformations in the literature that can handle singular (weak and strong) inte-
grands on a triangle domain. However, most of these transformations handle planar
triangles [19]. To integrate over curvilinear elements, as is the case in this thesis,
a first-order Taylor expansion of the surface is constructed around the singularity.
This tangential rectilinear domain may then be divided into triangular subdomains
and, then, integrated with quadrature rules constructed for planar triangles.

21



2.5.1 Subdivision of the quadrilateral element

The mapping r(u, v), given by (2.22), describes the curved element Se given the
constraints −1 ≤ u ≤ +1 and −1 ≤ v ≤ +1. Although it is possible to evaluate the
mapping outside the region −1 ≤ u ≤ +1 and −1 ≤ v ≤ +1, it should be avoided
since it involves extrapolation based on higher-order polynomials, where the result
often yields an extension of the curved surface that intersects with itself.

A method to alleviate this problem is to construct a first-order Taylor expansion
of the surface, where the linearized surface is the tangent plane evaluated at the
point r0, where r0 is the closest point on the surface to the fixed field point rf .
The linearization point r0 is determined by minimizing the cost function g(u, v) =
|r(u, v) − rf |2 subject to the constraints −1 ≤ u ≤ +1 and −1 ≤ v ≤ +1. By
choosing this cost function, the derivatives ∂g

∂u
and ∂g

∂v
are continuous everywhere

and a gradient-based solver may be exploited. The solution to the minimization
problem is (u0, v0) and we have r0 = r(u0, v0).

The linearization rlin(u, v) of the curved surface is constructed from the co-
variant basis vectors (a0, b0) evaluated at r0, which yields

rlin(u, v) = r0 + a0(u− u0) + b0(v − v0). (2.51)

In the following, we refer to rlin(u, v) as the the linearized surface that corresponds
to the curved surface r(u, v), where no constraints are applied to u and v. Also,
we introduce the the linearized element S̄e of the curved element Se, where the
linearized element S̄e is described by (2.51) subject to the constraints −1 ≤ u ≤ +1
and −1 ≤ v ≤ +1.

The linearized element S̄e is subdivided by constructing triangles based on the
corner vertics of S̄e and the projected point rp. The projection point rp is de-
termined by minimizing the cost function glin(u, v) = |rlin(u, v) − rf |2 without any
constraints on u and v. It is shown in Appendix B that the analytical solution to
the minimization problem is found by solving the system of linear equations |a0|2 a0 · b0

a0 · b0 |b0|2

 up

vp

 =

 a0 · (rf − r0 + a0u0 + b0v0)

b0 · (rf − r0 + a0u0 + b0v0)

 , (2.52)

where (up, vp) is the solution and we have rp = rlin(up, vp). Four possible cases of
the outcome of the subdivision are shown in Figures 2.6 - 2.9.

In practice, there are cases when rp is located very close to an edge or a corner
of S̄e such that extremely thin triangles are formed, where the height of the triangle
is much shorter than its base and the area of the triangle is close to zero. These
triangles are called slivers and the potential problem with slivers is discussed in later
chapters.

2.5.2 Mapping to suitable domain for product integration
rules

Thus, S̄e is sub-divided into a set of sub-triangles and, possibly, sub-quadrilaterals.
In the following, we consider a sub-triangle as a special case of a sub-quadrilateral.

22



-1

Y

0

1
-1

0

X

1

1.5

1

0.5

Z

(a) Curved surface Se and linearized surface
S̄e.

u

-1.5 -1 -0.5 0 0.5 1 1.5

v

-1

-0.5

0

0.5

1

(b) Corresponding sub-division of triangles
on the reference element.

Figure 2.6: Linearization point and projection point coincide inside the curved element. This
yields r0 = rp with −1 ≤ u ≤ +1 and −1 ≤ v ≤ +1. We get four sub-triangles on S̄e.

-0.5

Y

0

0.5-1

-0.5
X

0

0.5

1.5

1

0.5

Z

(a) Curved surface Se and linearized surface
S̄e.

u

-1.5 -1 -0.5 0 0.5 1 1.5

v

-1

-0.5

0

0.5

1

(b) Corresponding sub-division of triangles
on the reference element.

Figure 2.7: Linearization point and projection point coincide on an edge of the curved element.
This yields r0 = rp with |u| = 1 or |v| = 1, which corresponds to four special cases for the four
edges. We get three sub-triangles on S̄e.

We refer these sub-elements as generalized quadrilaterals, where the triangle is fea-
tured as a special case.

A number of numerical integration schemes in the existing literature are formu-
lated as a transformation from a triangle on S̄e to a quadrilateral domain in a new
coordinate system, which we refer to as the (p, q)-coordinate system below. It is
feasible to extend such a numerical integration procedure in order to apply to the
generalized quadrilaterals that we have given in the following observations:

– The projection point rp corresponds to p = q = 0. This implies that the
singularity of the Green’s function is most pronounced at the origin of the
(p, q)-coordinate system and the objective is to find a quadrature scheme that
efficiently deals with this singularity.

– The two (straight) edges of the generalized quadrilateral on S̄e that coincide

23



-1
-0.5

Y
0

0.5
-1

-0.5
X

0

0.5

1.5

1

0.5

Z

(a) Curved surface Se and linearized surface
S̄e.

u

-1.5 -1 -0.5 0 0.5 1 1.5

v

-1

-0.5

0

0.5

1

(b) Corresponding sub-division of triangles
on the reference element.

Figure 2.8: Linearization point and projection point coincide on an corner of the curved
element. This yields r0 = rp with |u| = |v| = 1, which corresponds to four special cases for the
four corners. We get two sub-triangles on S̄e.

-0.5

Y

0

0.5-1

-0.5
X

0

0.5

1.5

1

0.5

Z

(a) Curved surface Se and linearized surface
S̄e.

u

-1 -0.5 0 0.5 1 1.5

v

-1

-0.5

0

0.5

1

(b) Corresponding sub-division of triangles
on the reference element.

Figure 2.9: Linearization point on an edge or at a corner of the curved element and the
projection point is located elsewhere on the linearized surface. This yields r0 6= rp with |u| > 1
and/or |v| > 1 for rp, where the diagonals |u| = |v| are excluded as this case is treated in Fig. 2.8.
Notice that r0 is determined subject to the constraint −1 ≤ u ≤ +1 and −1 ≤ v ≤ +1 and, here,
at least one of these constraints is active. In addition, we have a number of possible cases for the
sub-division that may feature both triangles and quadrilaterals on S̄e.

with the two radial straight lines that intersect at rp are mapped to straight
edges in the (p, q)-plane, where these straight edges can be described by a
constant p-coordinate and q > 0. The two other straight edges of S̄e map onto
curves in the (p, q)-plane, where these curves may be expressed as functions
q = qmin(p) and q = qmax(p). If the coordinates (p, q) are chosen to be the
polar coordinates (φ, ρ), then a typical typical expression for these functions
could be ρ(φ) = h/ sinφ, where h is the distance from the projection point to
a line with constant u-coordinate in Figure 2.9(b) of the middle generalized
quadrilateral with four sides and φ is the polar angle.

24



– A generalized quadrilateral that specializes to a triangle on the linearized ele-
ment S̄e may take different forms in the (p, q)-plane. Here, the collapsed edge
of the generalized quadrilateral that corresponds to a point on S̄e is trans-
formed to a point in the (p, q)-plane.

Next, we construct quadrature schemes suitable for integration on the generalized
quadrilaterals expressed in the (p, q)-plane.

– The integration limits qmin(p) and qmax(p) are generally piecewise functions
with a discontinuous derivative at a small number of discrete points p1, p2, . . . , pN .

– We form a product rule for 2D by means of one-dimensional Gauss-quadrature
schemes, where all quadrature points fall on lines of constant p.

– First, we integrate with respect to p along segments pi < p < pi+1 and, then,
we integrate with respect to q from qmin(p) to qmax(p).

The quadrature weights and quadrature points are then mapped from the (p, q)-
plane, via the linearized element S̄e to the curved element Se.

This procedure is illustrated in the following. Consider an integral of the form

I =

∫
Se
h(ξ, ζ) dS, (2.53)

where ξ and ζ are coordinates on the curved element Se and h(ξ, ζ) represents
the product of a basis function and the scalar Green’s function. As described in
Section 2.2.4, the intregral over Se is transformed to the reference element according
to dS = D(u, v)dudv. However, the integration could also be performed over the
linearized element S̄e according to dSlin = D0dudv, where D0 = |a0×b0| is computed

at the point r0. This implies dS = D(u,v)
D0

dSlin.

The linearized surface element S̄e is then subdivided into generalized quadrilat-
erals and the integration over S̄e is converted to a sum over all resulting generalized
quadrilaterals Q, where the integration over each generalized quadrilateral is trans-
formed to the (p, q)-coordinate system according to dSlin = JTdqdp. Here, JT is the
Jacobian of the transformation rule, which cancels the singularity (weak or strong)
in the Green’s function. Hence, the integral (2.53) can finally be rewritten as∫

Se
h(ξ, ζ) dS =

∑
Q

∫ pmax

pmin

∫ qmax(p)

qmin(p)

h(u(p, q), v(p, q))
D(u, v)

D0

JTdqdp. (2.54)

2.5.3 Integrands with weak singularity

The term that by immidiate inspection possess a weak singularity (1/r) in the near-
field representation of the operators (1.43) and (1.44) is∫

S

G(r, r′)X(r′) dS ′. (2.55)

It is shown in [20] that the tangential component of the K-operator∫
S

([∇G(r, r′)]×X(r′))tan dS ′, (2.56)

25



also possess a weak singularity when the field point is located on the surface.
For integrands with a weak singularity, we exploit the Radial-Angular scheme

[21] and apply it according to the procedure outlined in Section 2.5.2. This approach
works for both the true and the near singularity.

2.5.4 Integrands with strong singularity

The terms that possess strong singularities (1/r2) in the near-field representation of
the operators (1.43) and (1.44) are∫

S

∇G(r, r′)∇′ ·X(r′) dS ′, (2.57)∫
S

[∇G(r, r′)]×X(r′) dS ′, (2.58)

where (2.58) contains a strong singularity for the near-singularity case.

True singularity

The integrand in (2.57) is the only term in the MoM formulation that possess a
true strong singularity. A method first developed by Guiggiani [22] for handling
true strong singularities in computational mechanics was reiterated by Weile [23]
for computational electromagnetics. In this scheme, the integral is divided into
parts that posses weak singularities and a part which is evaluated in the Cauchy
principal value sense. The complete derivation can be found in Appendix C, below
we just state the final results:

∫
S

∇G(r, r′)∇′ ·X(r′) dS ′ =

∫ 2π

φ=0

∫ β(φ)

ρ=0

[
Γ(ρ, φ)− γ(φ)

ρ

]
dρdφ

+

∫ 2π

φ=0

γ(φ) ln[β(φ)P (φ)]dρdφ (2.59)

+
∇ ·X(r)

2
n̂

The complete expressions for Γ, γ, P and β can be found in Appendix C.

Near singularity

For integrands with a near singularity, we exploit the scheme designed by Botha [11]
and apply it according to the procedure outlined in 2.5.2.

26



Chapter 3

Results

A numerical MoM-solver is implemented in MATLAB and it is tested in three differ-
ent settings: (i) scattered fields calculated from given source currents on a sphere;
(ii) induced currents and far-field solution for a PEC scatterer given an incident
plane wave; and (iii) induced currents and the near-field solution for two adjacent
PEC spheres illuminated by a plane wave.

3.1 Integration of known surface current

In this section, the scattered fields from a given source current are analyzed. This
serves as an intermediate test that verifies the expressions for the operators (1.37)
and (1.38). The analytical expressions for the surface currents are derived from the
multipole expansion, where a single mode is chosen for simplicity.

The MoM-formulation is tested by calculating the coefficients for each degree of
freedom of the current expansion (2.41) given the currents

J s(θ, φ) = θ̂

√
3

2π

E0

2Z0

1 + jka

(ka)2
e−jka sin(θ), (3.1)

M s(θ, φ) = φ̂

√
3

2π

E0

ka

[
1

2

(
kah

(2)
0 (ka)− h(2)

1 (ka)
)

sin(θ)

]
, (3.2)

where a is the radius of the sphere, E0 is a complex amplitude and h
(2)
1 , h

(2)
0 are

spherical Hankel functions of second kind. The surface currents (3.1) and (3.2)
produce scattered fields given by

E(r, θ, φ) =

√
3

2π

E0

kr

[
r̂h

(2)
1 (kr) cos(θ) + θ̂

1

2

(
h

(2)
1 (kr)− krh(2)

0 (kr)
)

sin(θ)

]
,

(3.3)

H(r, θ, φ) = −φ̂
√

3

2π

E0

2Z0

1 + jkr

(kr)2
e−jkr sin(θ). (3.4)

The expressions for the analytical surface currents and fields are derived in Appendix
D.

27



Given the surface currents, the scattered fields may be computed from the nu-
merical expressions for the current expanded in the known coefficients and basis
functions. The coefficients are computed from∫

S

wi ·Xa(r) dS =

∫
S

wi ·Xn(r) dS, (3.5)

where Xa is the analytical expression for the current and Xn is the numerical
expression for the current given by (2.41). Equation (3.5) is solved for the weighting
function (2.46) and the expression for the basis functions (2.41), which yields

XU
ij,e = (a ·X) det Je, (3.6)

XV
ij,e = (b ·X) det Je, (3.7)

where X is the analytical surface current J or M , a and b are the covariant basis
vectors.

3.1.1 Error analysis

In this section, the error of the MoM-formulation is analyzed, given a known source
current. The electric and magnetic field are computed by Eqs. (1.37), (1.38) and
compared to the analytical expressions (3.3) and (3.4). The absolute and the relative
error for a vector field F can be calculated according to

Absolute Error = |F num − F ana|, (3.8)

Relative Error =
|F num − F ana|
|F ana|

, (3.9)

where F num is the numerically computed result and F ana is the analytical field. We
asssume that total error etot can be expressed as

etot = emap + evec + eint, (3.10)

where emap is the error associated with the mapping of the geometry, evec is the error
associated with the discretization of the current and eint is the error that stems from
the integration scheme.

In the following tests, the error of the E-field and H-field are analyzed by
varying one of the spatial coordinates. The radius of the sphere is set to a = 1m,
the frequency is set by the ratio a/λ = 0.1, the order of the mapping and integration
scheme fixed to 12 and the order of the current discretization set to p = 2, 4, 6.
The sphere is discretized with six curvilinear quadrilateral cells. Thus, we expect
that evec is significantly larger than emap and eint.

For direct Gaussian quadrature, we compute the fields as a function of the field
point as it approaches the surface of the sphere. Here, the dyadic Green’s functions
(1.63) and (1.64) are used for the L and K operators. The relative and absolute
errors for the fields E and H are calculated with (3.9) and (3.8), where the results
are shown in Fig. 3.1.

28



r [m]

1.5 2 2.5 3

R
e
l
a
t
i
v
e
E
r
r
o
r
[
-
]

10
-4

10
-3

10
-2

10
-1

10
0 p = 2

p = 4

p = 6

(a) Relative error of E.

r [m]

1.5 2 2.5 3

A
b
s
o
l
u
t
e
E
r
r
o
r
[
V
/
m
]

10
-5

10
-4

10
-3

10
-2

10
-1

10
0

p = 2

p = 4

p = 6

norm(E)

(b) Absolute error of E.

r [m]

1.5 2 2.5 3

R
e
l
a
t
i
v
e
E
r
r
o
r
[
-
]

10
-4

10
-3

10
-2

10
-1

p = 2

p = 4

p = 6

(c) Relative error of H.

r [m]

1.5 2 2.5 3

A
b
s
o
l
u
t
e
E
r
r
o
r
[
A
/
m
]

10
-8

10
-7

10
-6

10
-5

10
-4

10
-3

p = 2

p = 4

p = 6

norm(H)

(d) Absolute error of H.

Figure 3.1: Errors of E and H by varying the distance of the field point in the direction
θ = 45◦, φ = 45◦; solid curve - p = 2; dashed curve - p = 4; and dashed-dotted curve - p = 6.

The errors increase as the field point moves towards the surface of the sphere.
For sufficiently small distance, direct Gaussian quadrature yields large errors due
to the 1/r dependence of the Green’s function. We also conclude that the error
decreases for higher order of p.

Next, the angular dependence of the error is analyzed, i.e. the field point is varied
at a fixed distance 3a around the sphere. The remaining paramters that describe
the problem are fixed to the same values that are used in Fig. 3.1. The relative
mean square angular error is computed by

Relative Error =

[
1

4π

∮
SΩ

|F num − F ana|2

|F ana|2
dΩ

] 1
2

, (3.11)

where dΩ is the solid angle and the fields are sampled at a distance r = 3a. The
integral (3.11) is evaluated with 12-th order Gaussian quadrature scheme in both θ-
and φ-direction. The result is shown in Fig. 3.2.

The graphs for E and H are nearly indistinguishable. A linear decrease of the
error in the above semi-log graph indicates an exponential convergence of the error
as p is increased.

Finally, the error is analyzed as a function of the frequency. The remaining

29



Basis function order p [-]
2 3 4 5 6

R
el
a
ti
v
e
E
rr
o
r
[-
]

10
-3

10
-2

10
-1 E-Field

H-Field

Figure 3.2: Plot showing the mean square angular error for E and H.

paramters that describe the problem are fixed to the same values that are used in
Fig. 3.1. However, the field point is fixed to a distance 3a from the origin. The
frequency is varied such that a/λ ∈ [10−4, 102] and the results are shown in Fig. 3.3.

For wavelengths λ > 100a, the problem tends to an electrostatic problem. Thus,
the error increases for the magnetic field for λ > 100a, whereas the error stays
constant for the electric field. For shorter wavelenghts, the error increases for both
the electric and magnetic field and the relative error is rather large for λ < a.

3.2 Scattering from PEC sphere

In this section, the MoM-formulation is tested for a PEC sphere of finite radius
a. The excitation source is assumed to be an incoming plane wave, where the
electric field given by Eq. (2.19). We set the propagation direction to k̂ = −ẑ and
the polarization to E0 = x̂E0 with unit amplitude E0 = 1. Given the incoming
field, the EFIE for the impedance matrix (2.13) and the excitation vector (2.15) are
computed. Next, the electric surface current is determined by solving the system
of linear equations (2.12). With the coefficients of the induced current known, the
scattered electric field is calculated with (1.37). In this case, the exact value for
the scattered electric far-field may be expressed as a Mie-series expansion [25] for
comparison. A quantity of interest is the radar cross section (RCS), which is defined
as [7]

σ = 4πr2 |E
s|2

|Ei|2
, (3.12)

where the scattered field Es is computed in the far-field region. The bistatic RCS
is computed in the φ = 0◦ plane and compared to the Mie-series expansion for a
frequency corresponding to a/λ = 0.1. The sphere is discretized with 24 curvilinear
quadrilateral cells. The RCS as a function of θ is plotted in Fig. 3.4. The EFIE

30



a/λ [-]

10
-4

10
-2

10
0

10
2

R
e
l
a
t
i
v
e
E
r
r
o
r
[
-
]

10
-4

10
-2

10
0

p = 2

p = 4

p = 6

(a) Relative error of E.

a/λ [-]

10
-4

10
-2

10
0

10
2

A
b
s
o
l
u
t
e
E
r
r
o
r
[
V
/
m
]

10
-8

10
-6

10
-4

10
-2

10
0

10
2

p = 2

p = 4

p = 6
norm(E

ana
)

(b) Absolute error of E.

a/λ [-]

10
-4

10
-2

10
0

10
2

R
e
l
a
t
i
v
e
E
r
r
o
r
[
-
]

10
-4

10
-2

10
0

10
2

p = 2

p = 4

p = 6

(c) Relative error of H.

a/λ [-]

10
-4

10
-2

10
0

10
2

A
b
s
o
l
u
t
e
E
r
r
o
r
[
V
/
m
]

10
-8

10
-6

10
-4

10
-2

10
0

p = 2

p = 4

p = 6
norm(H

ana
)

(d) Absolute error of Z0H.

Figure 3.3: Errors of E and H as a function of the normalized frequency. The distance of the
field point is fixed to r = 3a in the direction θ = 45◦, φ = 45◦; solid curve - p = 2; dashed curve -
p = 4; and dashed-dotted curve - p = 6.

result agrees very well with the Mie-series expansion, even for low-order vector basis
functions. The root-mean-square error (RMSE) of the RCS can be calculated with

RMSE =

[
1

Nθ

Nθ∑
i=1

|σEFIE(θi)− σMie(θi)|2

|σMie(θi)|2

] 1
2

, (3.13)

where Nθ is the number of point sampled in θ. The RMSE as function of vector basis
order is plotted in Fig. 3.5. Again, we note that the error decreases exponentially
with the order p of the divergence-conforming basis functions. Moreover, the bistatic
error is analyzed as a function of cell size h. The order of convergence in the cell size
hi (for i = 0, 1, 2) is estimated through carrying out computations for a geometric
sequence of cell sizes such that hi/hi+1 = hi+1/hi+2. We use cell sizes hi that
correspond to discretizing the surface of the sphere with 6, 24 and 96 curvelinear
quadrilateral elements. The bistatic error is calculated from Eq. (3.11) and the
result is shown in Fig. 3.6. For orders p = 2 and 4, the results form straight lines
in the log-log plot, with a steeper slope for the latter case. For the order p = 6, the
error increases unexpectedly from hi/h0 = 0.5 to hi/h0 = 0.25. The leading error

31



θ(deg)
0 50 100 150 200

R
C
S
[d
B
sm

]

-25

-20

-15

-10

-5

0

5

p = 2

p = 4

p = 6

Mie

Figure 3.4: RCS as a function of θ for a/λ = 0.1.

Basis function order p [-]

2 3 4 5 6

R
e
l
a
t
i
v
e
R
M
S
E

[
-
]

10
-6

10
-5

10
-4

10
-3

10
-2

10
-1

Figure 3.5: RMSE in the RCS of the EFIE solution for a/λ = 0.1.

can be assumed to be the lowest-order term in the series expansion of the error [8]

I(h) ≈ I0 + Iαh
α. (3.14)

where I0 is the extrapolated result and α is the order of convergence. For cell sizes
that form a geometric series, it is possible to estimate the order of convergence [8]
by

α = ln

[
I(hi)− I(hi+1)

I(hi+1)− I(hi+2)

]
/ ln

[
hi
hi+1

]
. (3.15)

32



hi/h0 [-]

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

R
e
l
a
t
i
v
e
E
r
r
o
r
[
-
]

10
-6

10
-4

10
-2

10
0

p = 2

p = 4

p = 6

Figure 3.6: Bistatic error of the EFIE solution as a function of the normalized cell size h for
a/λ = 0.1; solid curve - p = 2; dashed curve - p = 4; and dashed-dotted curve - p = 6.

When applied to the computed result in Fig. 3.6 for hi/h0 = 1, 0.5, and 0.25, this
formula gives α = 3.0318 for p = 2 and α = 4.6429 for p = 4. For p = 6 it is
not possible to get a realiable order of convergence from the numerical result. The
leading error is expected to scale as α = 2p. By inspection of the result for the cases
p = 2 and 4, it is noted that the resulting order of convergence is approximately a
factor of 0.75 and 0.56, respectively, of the theoretical value α = 2p.

3.3 Scattering from two PEC spheres

In this section, the MoM-formulation is tested on two adjacent PEC spheres, where
both sphere have the radius a. The centers of the spheres are separated by 2.5a
and they are illuminated by a plane wave with a frequency that corresponds to
a/λ = 0.1, where E0 = 1. The near-field around the spheres is calculated for (i)
horizontal polarization - an incoming electric field polarized along the straight line
between the centers of the spheres and (ii) vertical polarization - an incoming field
polarized perpendicular to the straight line between the centers of the spheres.

Figure 3.7 shows the electric field in the plane z = 0 for an incident electric
field with polarization along the straight line between the center of the sphere, i.e.
horizontal polarization as described above. In Fig. 3.7, it is possible to see an
enhanced field with roughly 2.5 times the amplitude of the incoming field between
the spheres.

Figure 3.8 shows the corresponding field plot for the incoming electric field ver-
tically polarized as described in (ii) above. The enhanced field between the spheres
does not occur for this polarization. Here, the maximum amplitude of the scattered
near-field is only a factor 1.6 times the amplitude of the incoming field, thus, it
is possible to conclude that the interaction between the spheres is stronger for the

33



Figure 3.7: The electric field around the spheres for a/λ = 0.1. The background color represents
the magnitude of the field and the red vector field the direction of the field.

Figure 3.8: The electric field around the spheres for a/λ = 0.1. The background color represents
the magnitude of the field and the red vector field the direction of the field.

horizontal polarization described in (i).

34



Chapter 4

Conclusions and future work

This thesis presents a Method of Moments (MoM) that exploits collocation in com-
bination with higher-order interpolatory divergence-conforming basis functions for
Maxwell’s equations. In this chapter, the main findings of the work are presented
together with some suggestions for future work.

4.1 Conclusions

The main focus is to investigate scattering problems treated by higer-order MoM.
The differential Maxwell’s equations are reformulated as the electric field integral
equations by introducing a Green’s function. The electric field integral equation
(EFIE) is reduced to a computational domian that only involves the surface of the
scatterer.

The MoM formulation is derived by the weighted residual method applied to
the EFIE. We exploit a collocation scheme for testing the EFIE, which could effec-
tively reduce the computational work required to evaluate the integrals in the EFIE
formulation.

The geometry is discretized by high-order curvelinear quadrilateral cells, where
Lagrangian polynomials are chosen for interpolation. An equal area mapping is used
to discretize the surface of the sphere into cells of uniform size.

The currents are expanded divergence-conforming basis functions, with inter-
polatory nodes that are collocated with the underlying quadrature schemes. This
reduces the integration for any given order to a single evaluation of the integrand
at the interpolatory node.

The singular integrals in the resulting EFIE formulation could effectively be
treated by subdividing the quadrilateral cell into triangles, where it is possible to
cancel the singularity with well-chosen coordinate transformations. However, trian-
gles resulting in slivers are hard to treat accurately with numerical integration and
they may produce errors that are large. Singularity extraction is used for strong
singularities.

The results show that the computer implemention of the MoM formulation could
effectively reproduce analytical results for PECs. The MoM formulation is showed
to converge exponentially with increasing vector basis order p. However, the current
computer implementation could not reproduce a leading error proportional to h2p,

35



where h is the cell size.

4.2 Future work

The topics covered in this work could be applied on many interesting problems. The
applicability of this work could be extended by including the modeling of dielectric
objects. For example, the field of plasmons is connected to internal oscillations of
electrons inside a metal. This corresponds to the presence of electric fields inside
the metal and the PEC approximation can not deal with those cases. A natural
continuation of this work would then be to include modeling of dielectrics.

36



Bibliography

[1] H. Jiang, L. Xu, and K. Zhan, “Joint tracking and classification based on
aerodynamic model and radar cross section,” Pattern Recognition, vol. 47, no. 9,
pp. 3096–3105, 2014.

[2] O. I. Sukharevsky, “Electromagnetic Wave Scattering by Aerial and Ground
Radar Objects,” no. 2, pp. 162–167, 2015.

[3] L. Ying, W. Zhe, H. Peilin, and L. Zhanhe, “A New Method for Analyzing
Integrated Stealth Ability of Penetration Aircraft,” Chinese Journal of Aero-
nautics, vol. 23, no. 2, pp. 187–193, 2010.

[4] Y. Li, J. Huang, S. Hong, Z. Wu, and Z. Liu, “A new assessment method for the
comprehensive stealth performance of penetration aircrafts,” Aerospace Science
and Technology, vol. 15, no. 7, pp. 511–518, 2011.

[5] H. Garcia, R. Sachan, and R. Kalyanaraman, “Optical Plasmon Properties of
Co-Ag Nanocomposites Within the Mean-Field Approximation,” Plasmonics,
vol. 7, no. 1, pp. 137–141, 2012.

[6] D. Albinsson, “Complex plasmonic nanostructures for materials science and
catalysis,” 2015. 64.

[7] J. D. Jackson, “Classical Electrodynamics,” 1999.

[8] A. Bondeson, T. Rylander, and P. Ingelström, Computational Electromagnetics.
External organization, 2005. 222.

[9] W. Gibson, The Method of Moments in Electromagnetics. No. 2, 2007.

[10] P. W. Fink, D. R. Wilton, and M. a. Khayat, “Simple and Efficient Numerical
Evaluation of Near-Hypersingular Integrals,” Antennas and Wireless Propaga-
tion Letters, IEEE, vol. 7, pp. 469–472, 2008.

[11] M. M. Botha, “Numerical Integration Scheme for the Near-Singular Green
Function Gradient on General Triangles,” IEEE Transactions on Antennas and
Propagation, vol. 63, no. 10, pp. 4435–4445, 2015.

[12] G. B. Arfken, Mathematical Methods for Physicists, 4th ed., vol. 67. 1999.

[13] R. F. Boisvert and D. W. Lozier, “Handbook of Mathematical Functions,” A
Century of Excellence in Measurements, Standards, and Technology - A Chron-
icle of Selected NBS/NIST Publications, 1901-2000, pp. 135–139, 2001.

37



[14] D. Roşca and G. Plonka, “Uniform spherical grids via equal area projection from
the cube to the sphere,” Journal of Computational and Applied Mathematics,
vol. 236, no. 6, pp. 1033–1041, 2011.

[15] A. F. Peterson, Mapped Vector Basis Functions for Electromagnetic Integral
Equations, vol. 1. 2006.

[16] F. B. Hildebrand, Introduction to numerical analysis. 1987.

[17] A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics. 2007.

[18] S. Rao, D. Wilton, and a. Glisson, “Electromagnetic scattering by surfaces of
arbitrary shape,” IEEE Transactions on Antennas and Propagation, vol. 30,
no. 3, pp. 409–418, 1982.

[19] D. R. Wilton, F. Vipiana, and W. a. Johnson, “Evaluating Singular, Near-
Singular, and Non-Singular Integrals on Curvilinear Elements,” Electromagnet-
ics, vol. 34, no. 3-4, pp. 307–327, 2014.

[20] S. D. Gedney, “On Deriving a Locally Corrected Nyström Scheme from a
Quadrature Sampled Moment Method,” IEEE Transactions on Antennas and
Propagation, vol. 51, no. 9, pp. 2402–2412, 2003.

[21] M. a. Khayat and D. R. Wilton, “Numerical evaluation of singular and near-
singular potential integrals,” IEEE Transactions on Antennas and Propagation,
vol. 53, no. 10, pp. 3180–3190, 2005.

[22] M. Guiggiani and a. Gigante, “A General Algorithm for Multidimensional
Cauchy Principal Value Integrals in the Boundary Element Method,” Journal
of Applied Mechanics, vol. 57, no. 4, p. 906, 1990.

[23] D. S. Weile and X. Wang, “Strong singularity reduction for curved patches
for the integral equations of electromagnetics,” IEEE Antennas and Wireless
Propagation Letters, vol. 8, no. 2, pp. 1370–1373, 2009.

[24] Michael G. Duffy, “Quadrature Over a Pyramid or Cube of Integrands with a
Singularity at a Vertex,” SIAM Journal on Numerical Analysis, vol. 19, no. 6,
pp. 1260–1262, 1982.

[25] C. a. Balanis, Advanced Engineering Electromagnetics, vol. 52. 1989.

38



Appendix A

Far-field derivation for L and K

In this section, we want to find the far-field expressions for the operators

(LX)(r) =

[
1 +

1

k2
∇∇ ·

] ∫
S

G(r, r′)X(r′) dS ′, (A.1)

(KX)(r) = ∇×
∫
S

G(r, r′)X(r′) dS ′s, (A.2)

given the 3D, homogeneous, scalar Green’s function

G(r, r′) =
e−jk|r−r

′|

4π|r − r′|
. (A.3)

We start with the expression for L. We rewrite the distance |r−r′| between the
source point r′ and observation point r as a scalar product

|r − r′| =
√

(r − r′) · (r − r′) =
√
r2 − r′2 − 2r · r′

where r = |r| and r′ = |r′|. We can approximate the leading contribution from this
distance as

|r − r′| = r

√
1 + (

r′

r
)2 − 2r̂ · r

′

r
= r

{
1 +

1

2

[(
r′

r

)2

− 2r̂ · r
′

r

]
+ . . .

}
= r − r̂ · r′ +O(d2/r), as→∞,

(A.4)

where we have used
√

1 + x = 1 + x/2 + . . . and defined

d = max|r′|

as the maximum extension of the scatterer. The Green’s function may then be
rewritten as

G(r, r′) =
e−jk|r−r

′|

4π|r − r′|
=

exp(−jk(r − r̂ · r′ +O(d2/r))

4πr(1 +O(d/r))

=
e−jkr

4πr
e−jkr̂ ·r

′
(1 +O(kd2/r))(1 +O(d/r)).

39



The leading contribution from (A.1) becomes then

(LX)(r) =

[
1 +

1

k2
∇∇ ·

]
e−jkr

kr

k

4π

∫
S

e−jkr̂ ·r
′
X(r′) dS ′, (A.5)

where we have used r � d and r � kd2. It is convenient to define

P (r̂) =
k

4π

∫
S

e−jkr̂ ·k̂X(r′) dS ′, (A.6)

where P is essentially the spatial Fourier transform of the current density X eval-
uated in the point kr̂. Equation (A.5) can then be rewritten as

(LX)(r) =

[
1 +

1

k2
∇∇ ·

](
e−jkr

kr
P (r̂)

)
.

Evaluating the contribution from divergence operator gives

∇ ·
(
e−jkr

kr
P (r̂)

)
=
e−jkr

kr
∇ · P (r̂) + P (r̂) · ∇

(
e−jkr

kr

)
.

Since P (r̂) = P (φ, θ), we get in spherical coordinates

∇ · P (r̂) =
1

r sin θ

∂(sin θPθ)

∂θ
+

1

r sin θ

∂Pφ
∂φ
∝ 1

r
,

∇
(
e−jkr

kr

)
= r̂(−jk − 1

r
)
e−jkr

kr
.

This gives

1

k
∇ ·

[
e−jkr

kr
P (r̂)

]
=
e−jkr

(kr)2

1

sin θ

[
∂(sin θPθ)

∂θ
+
∂Pφ
∂φ

]
+ P (r̂) · r̂(−jk − 1

r
)
e−jkr

kr

= −jP (r̂)
e−jkr

kr
(1 +O((kr)−1)).

Similary, we get

1

k
∇
[

1

k
∇ ·

(
e−jkr

kr
P (r̂)

)]
=

1

k
r̂
∂

∂r

[
−jr̂ · P (r̂)

e−jkr

kr
(1 +O((kr)−1))

]
=

1

k
r̂

[
−jr̂ · P (r̂)(−jk − 1

r
)
e−jkr

kr
(1 +O((kr)−1))

]
= −r̂

[
r̂ · P (r̂)

e−jkr

kr
(1 +O((kr)−1))

]
Hence, the operator L may be rewritten as

(LX)(r) = [P (r̂)− r̂(r̂ · P (r̂))]
e−jkr

kr

= −r̂ × (r̂ × P (r̂))
e−jkr

kr
,

(A.7)

40



where we have used r � λ.
Following the above example for the operator K, we have

(KX)(r) = ∇×
∫
S

G(r, r′)X(r′) dS ′s

=
1

k
∇×

[
e−jkr

kr
P (r̂))

]
.

Evaluating the cross product in spherical coordinates gives

1

k
∇×

[
e−jkr

kr
P (r̂)

]
=

1

k

[
e−jkr

kr
∇× P (r̂) +∇

(
e−jkr

kr

)
× P (r̂)

]
=
e−jkr

(kr)2

[
r̂

1

sin θ

(
∂(Pφ sin θ)

∂θ
− ∂Pθ

∂φ

)
+ θ̂

1

sin θ

∂Pr
∂φ
− φ̂∂Pr

∂θ

]
+

[
r̂

1

k

(
−jk − 1

r

)
e−jkr

kr

]
× P (r̂)

= −jr̂ × P (r̂)
e−jkr

kr
(1 +O((kr)−1)).

Using r � λ, we finally arrive at

(KX)(r) = −jr̂ × P (r̂)
e−jkr

kr
. (A.8)

It is clear from (A.7) and (A.8) that the dominating term scales as 1/r and that the
far-fields possess no components in the direction of propagation r̂ outward from the
scatterer.

41



42



Appendix B

Projection point on tangent plane

Assume we have a curved surface described by the mapping r(u, v). A first-order
Taylor approximation of the curved surface around the point (u0, v0) is described by

rlin(u, v) = r0 + a0(u− u0) + b0(v − v0), (B.1)

where r0 = r(u0, v0) and a0 = ∂r
∂u

∣∣∣
(u0,v0)

, b0 = ∂r
∂v

∣∣∣
(u0,v0)

are the co-variant basis

vectors. Given a known field point rf , the projected point rp = r(up, vp) of rf onto
the linearized surface rlin(u, v) is found by minimizing the cost function

g(u, v) = |rlin(u, v)− rf |2 = |rlin(u, v)|2 − 2rlin(u, v)rf + |rf |2. (B.2)

The analytical solution is found by setting the partial derivatives to zero
∂g
∂u

∣∣∣
(up,vp)

= −2rf · a0 + 2rlin(up, vp) · a0 = 2a0 · (rlin(up, vp)− rf) = 0

∂g
∂v

∣∣∣
(up,vp)

= −2rf · b0 + 2rlin(up, vp) · b0 = 2b0 · (rlin(up, vp)− rf) = 0
, (B.3)

{
a0 · (r0 + a0(up − u0) + b0(vp − v0)− rf) = 0

b0 · (r0 + a0(up − u0) + b0(vfflp − v0)− rf) = 0
. (B.4)

Inserting the expression for rlin and collecting terms we arrive at the system of linear
equations  |a0|2 a0 · b0

a0 · b0 |b0|2

 up

vp

 =

 a0 · (rf − r0 + a0u0 + b0v0)

b0 · (rf − r0 + a0u0 + b0v0)

 , (B.5)

where (up, vp) is the solution.

43



44



Appendix C

Strong true singularity treatment

Consider the integral

I(r) =

∫
S

∇G(r, r′)∇′ ·X(r′) dS ′, (C.1)

where we have a true strong singularity located on the surface S. In the following,
we denote ∇′ ·X(r′) = s(r′). Then

I(r) = −
∫
S

∇
(
e−jk|r−r

′|

4π|r − r′|

)
s(r′) dS ′

=

∫
S

[
(1 + jk|r − r′|e−jk|r−r′|

4π|r − r′|2
r − r′

|r − r′|

]
s(r′) dS ′

=

R = |r − r′|

R̂ = (r − r′)/|r − r′|

=

∫
S1

[
(1 + jkRe−jkR

4πR2
R̂

]
s(r′) dS ′

+

∫
S2

[
(1 + jkRe−jkR

4πR2
R̂

]
s(r′) dS ′

= I1(r) + I2(r)

where S2 is a disk with the center located at the singulary and with radius ε, and
S1 = S ∩ S2. Assume that ε is so small that S2 is flat and perform integration in
cylindrical coordinates. Also, assume that ε is sufficiently small to assume statics,

45



i.e kR→ 0 so that 1 + jkR→ 1 and e−jkR → 1.

I2(r) '
∫
S2

[
1

4πR2
R̂

]
s(r′) dS ′

= {ρ(r′) ' constant on S2 ⇒ s(r)}

' ρ(r)

4π

∫ 2π

φ′=0

∫ ε

r′=0

[
−r̂(φ′)r′ + ẑz

[(r′)2 + z2]3/2

]
r′dr′dφ′

=
ρ(r)

2

∫ ε

r′=0

ẑzr′

[(r′)2 + z2]3/2
dr′

=
ρ(r)

2
ẑz

(
1

z
− 1√

ε2 + z2

)
lim
ε→0

I2(r) =
s(r)

2
ẑ =

s(r)

2
n̂

We are left with

I(r) =

∫
S1

[
(1 + jkR)e−jkR

4πR2
R̂

]
s(r′) dS ′ +

s(r)

2
n̂ (C.2)

For S1, we use the mapping r(u, v) described in (2.22) with corresponding covariant
basis vectors a, b (2.28) and D(u, v) = |a × b|. Let (u0, v0) be centered at S2 and
use the polar coordinates defined by

u = u0 + ρ cosφ

v = v0 + ρ sinφ

The excluded area S2 is bounded by |r(u, v)− r(u0, v0)| = ε. Since ε is very small,
we can use Taylor expansion

r(u, v) = r(u0, v0)+
∂r

∂u

∣∣∣
(u0,v0)

(u− u0) +
∂r

∂v

∣∣∣
(u0,v0)

(v − v0) . . .

⇒ |r(u, v)− r(u0, v0)| ' |a0(u− u0) + b0(v − v0)|
= |a0ρ cosφ+ b0ρ sinφ|
= ρ|a0 cosφ+ b0 sinφ| = ε.

Thus, we have S2 described by

ρε ≤
ε

|a0 cosφ+ b0 sinφ|
. (C.3)

The outer boundary is described in terms of straight line segments form a point
(u1, v1) to another point (u2, v2), which gives at the boundary

β(ξ) = ξ(u1, v1) + (1− ξ)(u2, v2) for 0 ≤ ξ ≤ 1.{
p(ξ) = ξu1 + (1− ξ)u2 − u0 = ρ cosφ

q(ξ) = ξv1 + (1− ξ)v2 − v0 = ρ sinφ

46



{
ρβ(ξ) =

√
p2(ξ) + q2(ξ)

φ(ξ) = arctan(p(ξ)/
√
p2(ξ) + q2(ξ))

or 
ρβ(φ) = u0(v1−v2)+u1(v2−v0)+u2(v0−v1)

(v2−v1) cosφ+(u1−u2) sinφ

ξ(φ) = (v2−v0) cosφ+(u0−u2) sinφ
(v2−v1) cosφ+(u1−u2) sinφ

,

which is useful in the following. We now have

I1(r) =

∫
S1

[
(1 + jkR)e−jkR

4πR2
R̂

]
s(r′) dS ′

=

∫ 2π

φ′=0

∫ ρ′β(φ′)

ρ′ε(φ
′)

[
(1 + jkR)e−jkR

4πR2
R̂

]
s(r′)D(u′, v′)ρ′ dρ′dφ′

=

{
Γ(u′, v′) = Γ(ρ′, φ′) =

[
(1 + jkR)e−jkR

4πR2
R̂

]
s(r′)D(u′, v′)ρ′

}
=

∫ 2π

φ′=0

∫ ρ′β(φ′)

ρ′ε(φ
′)

Γ(ρ′, φ′) dρ′dφ′

where 
r = r(u0, v0)

r′ = r(u′, v′)

R = |r(u0, v0)− r(u′, v′)|
R̂ = (r(u0, v0)− r(u′, v′))/|r(u0, v0)− r(u′, v′)|.

The integrand Γ exhibits a weak singularity (∝ 1/ρ′) as the point (u0, v0) is ap-
proached. The proportionality constant for the part of the integrand that scales as
1/ρ is given by

γ(φ′) = lim
ρ′→0

ρ′Γ(ρ′, φ′)

=


R̂ = r(u0,v0)−r(u′,v′)

|r(u0,v0)−r(u′,v′)| = Û0(u0−u′)+V̂ 0(v0−v′)
|Û0(u0−u′)+V̂ 0(v0−v′)|

R = |r(u0, v0)− r(u′, v′)| = |Û 0(u0 − u′) + V̂ 0(v0 − v′)|
Static approx:(1 + jkR)→ 1 and e−jkR → 1

= lim
ρ′→0

−(a0(u0 − u′) + b0(v0 − v′))
4π|a0(u0 − u′) + b00(v0 − v′)|3

s(r′)D(u′, v′)(ρ′)2

= lim
ρ′→0

−(ρ′(a0 cosφ′ + b00 sinφ′))

4π(ρ′)3|a0 cosφ′ + b0 sinφ′|3
s(r′)D(u′, v′)(ρ′)2

=
1

4π

−(a0 cosφ′ + b0 sinφ′)

|a0 cosφ′ + b0 sinφ′|3
s(r′(u0, v0))D(u0, v0).

47



Thus, we have

I1(r) =

∫ 2π

φ′=0

∫ ρ′β(φ′)

ρ′ε(φ
′)

(
Γ(ρ′, φ′)− γ(φ′)

ρ′

)
dρ′dφ′

+

∫ 2π

φ′=0

∫ ρ′β(φ′)

ρ′ε(φ
′)

γ(φ′)

ρ′
dρ′dφ′

= Ireg
1 (r) + Ising

1 (r),

where

Ising
1 (r) =

∫ 2π

φ′=0

γ(φ′) ln

(
ρ′β(φ′)

ρ′ε(φ
′)

)
dφ′

=

∫ 2π

φ′=0

γ(φ′) ln

(
ρ′β(φ′)

|a0 cosφ+ b0 sinφ|
ε

)
dφ′

=

∫ 2π

φ′=0

γ(φ′) ln
(
ρ′β(φ′)|a0 cosφ+ b0 sinφ|

)
dφ′ − ln ε

∫ 2π

φ′=0

γ(φ′) dφ′

=

∫ 2π

φ′=0

γ(φ′) ln
(
ρ′β(φ′)|a0 cosφ+ b0 sinφ|

)
dφ′,

where the last integral vanished since γ(r′) is periodic on the interval [0, 2π] with a
zero mean. Finally, we get (when ε→ 0+)

I(r) =

∫ 2π

φ′=0

∫ ρ′β(φ′)

ρ′ε(φ
′)→0+

(
Γ(ρ′, φ′)− γ(φ′)

ρ′

)
dρ′dφ′

+

∫ 2π

φ′=0

γ(φ′) ln
(
ρ′β(φ′)|a0 cosφ+ b0 sinφ|

)
dφ′

+
s(r)

2
.

The first term in the above expression is regular because γ is defined as: γ(φ) =
limρ→0 ρΓ(ρ, φ), and may be integrated with direct Gaussian quadrature. The second
term is a regular contour integral and the last term is the Cauchy principal value
term.

48



Appendix D

Analytic solutions for a test
problem

Consider a sphere of radius a with surface currents J and M . For r > a, the
electromagnetic field from one mode in the multipole expansion are given in spherical
coordinates as

E(r, θ, φ) =

√
3

2π

E0

kr

[
r̂h

(2)
1 (kr) cos(θ) + θ̂

1

2

(
h

(2)
1 (kr)− krh(2)

0 (kr)
)

sin(θ)

]
,

(D.1)

H(r, θ, φ) = −φ̂
√

3

2π

E0

2Z0

1 + jkr

(kr)2
e−jkr sin(θ), (D.2)

where h
(2)
1 , h

(2)
0 are spherical Hankel functions of second kind. Further, we let E = 0

and H = 0 for r < a. Thus, we may calculate the currents J and M on the surface
r = a from the boundary conditions,

−n̂ ×E1 = M s, (D.3)

n̂ ×H1 = J s, (D.4)

n̂ ·D1 = ρs, (D.5)

n̂ ·B1 = %s. (D.6)

On the surface r = a of the sphere, we have n̂ = r̂ and Eqs. (D.3) and (D.4) yield

J s(θ, φ) = θ̂

√
3

2π

E0

2Z0

1 + jka

(ka)2
e−jka sin(θ), (D.7)

M s(θ, φ) = φ̂

√
3

2π

E0

ka

[
1

2

(
kah

(2)
0 (ka)− h(2)

1 (ka)
)

sin(θ)

]
, (D.8)

where we have used the expressions for the scattered fields (D.1), (D.2) and the
relations r̂ × r̂ = 0, r̂ × θ̂ = φ̂.

49


	Introduction
	Maxwell's equations
	Boundary conditions at interfaces between different media

	Integral representations of the electromagnetic fields
	Alternative representations (suitable for the near-field region)
	Far-field representations
	Green's functions


	Method of Moments
	Integral equations for PEC scatterers
	Representation of the geometry
	Discretization of a sphere
	Reference element and higher-order mapping
	Co- and contravariant basis vectors
	Surface integration

	Representation of the current density
	Divergence conforming basis functions

	Weighting functions – Petrov-Galerkin's method
	Singular integrals
	Subdivision of the quadrilateral element
	Mapping to suitable domain for product integration rules
	Integrands with weak singularity
	Integrands with strong singularity


	Results
	Integration of known surface current
	Error analysis

	Scattering from PEC sphere
	Scattering from two PEC spheres

	Conclusions and future work
	Conclusions
	Future work

	Bibliography
	Far-field derivation for L and K
	Projection point on tangent plane
	Strong true singularity treatment
	Analytic solutions for a test problem