# Development of a symplectic scheme with optimized numerical dispersion-relation equation to solve Maxwell's equations in dispersive media.

1. INTRODUCTIONWith the advent of ever-improving hardware in CPU/GPU computational platforms and numerical methods in parallel environment, computational electromagnetics has played an essential role in the design of modern optical devices applied in different areas. This computing alternative enables also an exploration of propagation details in the optical media of practical relevance. We are therefore motivated to develop a physically correct and numerically efficient solver for more accurately predicting the electromagnetic (EM) wave propagation in the frequency-dependent and frequency-independent optical media.

Prediction of EM wave propagation in optical media can be carried out either in frequency or in time domain. In this study the finite difference method is adopted to solve the Maxwell's equations in time domain so that the evolution of wave can be more clearly understood. Within the FDTD context, approximation of the time and spatial derivative terms in the Ampere's and Faraday's equations must be physically correct and numerically efficient. For the Maxwell's equations defined in vacuum, for example, the geometric symplectic structure must be retained all the time. Otherwise, the Hamiltonian and Casimir in this Hamiltonian differential system can no longer be preserved in particular after a long time simulation [1]. Violation of the symplecticity-preserving property prohibits properly designing many long-range wave propagation optical devices. Care on the approximation of first-order spatial derivative terms must be also taken into a serious consideration to avoid the dispersion error. Such an error type can result in a complete erroneous propagation speed, thereby leading to an unphysically oscillatory solution. Besides the errors generated from the approximation of the time and space derivative terms, calculation of the Ampere's and Faraday's equations is subject to the satisfaction of discrete Gauss's law [2]. The above three theoretical considerations prompted our previous development of the dispersion-relation, divergence-free, and symplecticity triple-preserving scheme for solving the Maxwell's in a computationally more simple and a programmingly more easy non-staggered (or colocated) grids [3].

In this study, we attempt to extend our previously developed scheme for solving the Maxwell's equations in free space to simulate the frequency-dependent media. Three dispersive medium types, which comprise the Debye, Lorentz and Drude dispersive media, are all included in the generalized code for the Maxwell's equations. To numerically bridge the time increment and the spatial interval in the TM wave prediction, the means of minimizing the difference between the dispersion relation equation and the exact dispersion relation equation is adopted in this study rather than the minimization of modified wave number derived in our previous paper [3,4]. The necessity of deriving the numerical dispersion relation equation in the current minimization of dispersive error prompts us to adopt the explicit partitioned Runge-Kutta non-iterative scheme rather than the implicit Runge-Kutta symplectic scheme proposed in [4]. In addition, to satisfy the discrete Gauss law the Yee's staggered grid approach [5] is adopted rather than the non-staggered grid approach published earlier in [3, 4] to simplify the derivation of numerical dispersion relation equation.

To avoid an erroneous wave reflection from the truncated boundary, one can either prescribe a proper radiation boundary condition along the truncated boundary itself or attach a finite width layer at a location that is immediately adjacent to the physical boundary to absorb the unphysical wave reflection. In this study the CPML (Convolution Perfect Matching Layer) of Roden and Gedney [6] rather than the PML of Berenger [7] is adopted since the former matching layer approach has been reported to be able to absorb the wave more effectively.

The rest of the paper is organized as follows. In Section 2, the working equations applicable to non-dispersive and dispersive media are developed in free space as well as in perfect matching layer. In Section 3, we will present a splitting solution algorithm so that the core of the Maxwell's equations can be rigorously approximated not only in space but also in time. In Section 4, the explicit partitioned Runge-Kutta symplectic temporal scheme, which involves no iterative procedure, is proposed to preserve the Hamiltonian in the ideal Maxwell's equations. This explicit scheme is also essential for us to derive the numerical dispersion relation equations in Section 5 for the cases investigated in one and two dimensional domains. In Section 6, we analytically validate the proposed numerical methods and then in Section 7 to address the distinction between the three investigated dispersive media. Finally, some concluding remarks are drawn in Section 8.

2. WORKING EQUATIONS

Most metals and many other materials exhibit their own dispersive characteristics. The permittivity and permeability of these dispersive media are functions of the optical frequency. For simplicity, in the current study only the electric permittivity is assumed to be frequency-dependent. Extension of the simulation to the case of frequency-dependent magnetic permeability and electric permittivity will be our future study.

Numerical simulation of electromagnetic fields will be carried out in time domain. While most constitutive equations for the dispersive media are available in frequency domain, we should adopt its time-domain counterpart so as to be consistent with the FDTD simulation. As a result, the electric permittivity [epsilon](t) becomes the sum of the relative permittivity [[epsilon].sub.r]([equivalent to] [[epsilon].sub.[infinity]][delta](t) + [chi](t)) considered at the infinite frequency [[epsilon].sub.[infinity]] and the electric susceptibility [chi](t) [8]

[epsilon](t) = [[epsilon].sub.o] ([[epsilon].sub.[infinity]][delta](t) + [chi](t)). (1)

In the above, t denotes time, [[epsilon].sub.0] is the permittivity in free space and [delta](t) stands for the delta function.

For a dispersive medium whose magnetic permeability is frequency independent, the Ampere's and Faraday's laws can be respectively represented below in time domain for modeling the time-varying electric field variable [E.bar] and the magnetic field variable [H.bar]

[[partial derivative]/[partial derivative]t]([epsilon]([x.bar],t) * [E.bar]([x.bar],t)) = [nabla] x [H.bar] - [[J.bar].sub.d], (2)

[mu][[[partial derivative][H.bar]]/[partial derivative]t] = -[nabla] x [E.bar]. (3)

In Equation (2), [epsilon] is equal to [[epsilon].sub.0][[epsilon].sub.r] while [mu] in Equation (3) is identical to [[mu].sub.0][[mu].sub.r]. The optical properties [[epsilon].sub.0], [[epsilon].sub.r], [[mu].sub.0] and [[mu].sub.r] represent the free space electric permittivity, relative electric permittivity, free-space magnetic permeability, and relative magnetic permeability, respectively. [[J.bar].sub.d] in Equation (2) denotes the polarization current and depends on the types of the investigated dispersive optical media. For simplicity, both of the volume electric and magnetic current densities are assumed to be zero under the source-free condition. The resulting differential system will be solved subject to the Gauss's law which comprises [nabla] x [B.bar] = [nabla] x [D.bar] = 0. The electric permittivity [epsilon] in (1) is a function of frequency for the sake of simplicity when performing the current simulation. The notation "*" in Equation (2) denotes the convolution operator and it is defined as f(t) * g(t) = [[integral].sub.0.sup.t]f(t - [tau])g([tau]) d[tau] for any two given functions f (t) and g(t).

In the literature three constitutive models are often applied to perform the FDTD electromagnetic wave simulations in Debye, Lorentz and Drude media. For these three distinct media, their susceptibility functions [chi]([omega]) vary with time t and frequency [omega]. As a result, the relative electric permittivity for the optical media having P number of poles can be expressed in terms of the frequency and they are summarized in Table 1. In this table, u(t) denotes the unit step function.

The parameters of the three dispersive media are chosen to be [[epsilon].sub.[infinity]] = 7, [[epsilon].sub.s] = 10, [tau] = 7.0 x [10.sup.-10]s for the Debye medium, [[epsilon].sub.[infinity]] = 1.5, [[epsilon].sub.s] = 3, [[omega].sub.p] = 40[pi]GHz, [[delta].sub.p] = 4[pi] GHz for the Lorentz medium and [[epsilon].sub.[infinity]] = 1, [f.sub.p] = 28.7 GHz, [[omega].sub.p] = 2[pi][f.sub.p], [[gamma].sub.p] = 20 GHz for the Drude medium [9]. The equations used for modeling the polarization current in the three investigated dispersive media are given below

(1) Debye medium

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

(2) Lorentz medium

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

(3) Drude medium

[[gamma].sub.p][[partial derivative][[J.bar].sub.d]/[partial derivative]t] + [[[partial derivative].sup.2][[[J.bar].sub.d][partial derivative][t.sup.2]] = [[epsilon].sub.0][[omega].sup.2.sub.p][[partial derivative][E.bar]/[partial derivative]t. (6)

In the above, [DELTA][epsilon] = [[epsilon].sub.s] - [[epsilon].sub.[infinity]], where [epsilon]s is the static or zero-frequency relative permittivity, [[epsilon].sub.[infinity]] the permittivity at infinite frequency, [tau] the pole relaxation time, [omega] the angular frequency, [[omega].sub.p] the undamped resonant frequency of the medium, and [[gamma].sub.p] the inverse of the pole relaxation time.

When simulating the electromagnetic wave propagation, it is normally necessary to truncate the analysis domain somewhere in free space so as to make the analysis computationally feasible. Truncation of physical domain leads very often to a wave reflection back to the domain of interest. Such an unphysical wave reflection will possibly interact with the incident physical wave. To circumvent this problem, one can either prescribe a proper set of differential equations on the truncated boundary or attach a layer that is capable of absorbing the reflected wave. In this study we adopt the means of attaching a perfect matching layer at a location that is immediately downstream of the truncated boundary.

Following the work of Berenger [7], the Maxwell's equations in (2)(3) are first rewritten into their equivalent time-harmonic equations. Besides the introduction of stretched coordinate metrics into the Maxwell's equations, the convolution operator is also accelerated using the recursive convolution (RC) method proposed originally by Lubber et al. [10]. The resulting modified Maxwell's equations, which take into account wave absorption, in frequency domain are then transformed back to the time domain equations by virtue of the Fourier transform. The Maxwell's equations in the convolution perfect matching layer (CPML) can then be derived.

In the calculation of the resulting electromagnetic equations in time domain, the computational challenge lies in the implementation of convolution terms shown in the last two terms on the right hand sides

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

Note that [k.sub.i] (i = x, y, z) is the wave number along the i direction and [[zeta].sub.w](t) * [partial derivative][H.sub.v](t)/[partial derivative]w or [[zeta].sub.w](t) * [partial derivative][E.sub.v](t)/[partial derivative]w (w = x,y,z; v = x,y,z) are the convolution terms, which are present only in PML regions.

By applying a piecewise constant approximation on the convolution terms shown in (7)-(8), the Ampere's and Faraday's equations in the convolution perfect matching absorbing layers are as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

In the above, [psi][E.sub.w,v] or [psi][H.sub.w,v] is defined as [[zeta].sub.w](t) * [partial derivative][E.sub.v](t)/[partial derivative]w or [[zeta].sub.w](t) * [partial derivative][H.sub.v](t)/[partial derivative]w. The TM-mode Maxwell's equations for a lossless dispersive medium under current investigation are simplified as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Note that the coefficient [b.sub.w] and [c.sub.w] (w = x or y) are nonzero only in PML region with normal interface boundaries. In the above, [b.sub.w] and [c.sub.w] (w = x or y) are both chosen to be the exponential functions given below

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Note that the value of w falls between 0 and d. In PML, [[sigma].sub.w] is increased from 0 at w = 0 to [[sigma].sub.w,max] at w = d. Similarly, in PML the value of [k.sub.w] is increased from 1 at w = 0 to [k.sub.w,max] at w = d [8].

3. SOLUTION ALGORITHM

In Section 2 we have derived the governing equations for the field variables [E.sub.z], [H.sub.x] and [H.sub.y] in a domain of two dimensions. This set of TM-mode equations applicable to the three investigated dispersive media, which are characterized by their constitutive equations for the relative electrical permittivities [[epsilon].sub.r], contains the ideal Maxwell's equations and the terms due to the polarization current [[J.bar].sub.d] summarized in Equations (4)-(6), and the absorbing terms defined only in the attached convolution perfect matching layer. The quality of simulating the EM wave propagation in a dispersive medium depends therefore on the prediction quality of the Maxwell's equations derived in free space, the suitability of the employed constitutive equations, which are tabulated in Table 1, the simulation quality on the Equations (4)-(6) for the polarization current, and the ability of the convolution terms added to absorb the possibly reflected waves from the truncated boundary.

In the above light, we are led to know that the entire set of equations in (11) can be decomposed into the components accounting for the propagation kernel, which denotes the ideal Maxwell's equations, the constitutive laws for the dispersive media, the polarization current of the optical media, and the wave absorption. The possibility of decomposing the equations into the specific components motivates us to numerically treat them separately using their respective suitable numerical methods. Following this line of thought, we propose in this paper the numerical methods to calculate the absorption terms [psi][H.sub.x,y], [psi][H.sub.y,x], [psi][E.sub.z,x], and [psi][E.sub.z,y] in Equations (4)-(6) and the polarization currents. In addition to these two calculations, the solution quality for the wave simulation in dispersive media depends highly on the way of discretizing the following Maxwell's equations in vacuum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)

We will address in next section that the symplectic structure in the equations can be retained in the course of approximating the equations in (13). Besides choosing the temporal discretization to preserve the symplectic nature in the Hamiltonian differential system, we will also optimize the dispersion relation equation when discretizing the first-order spatial derivative terms. It is best hoped that the chosen values of [DELTA]t and [DELTA][x.bar] are well paired so as to get the best numerical dispersion relation equations for the Ampere's and Faraday's equations.

4. EXPLICIT PARTITIONED RUNGE-KUTTA TEMPORAL DISCRETIZATION

Faraday's and Ampere's equations in an ideal medium constitute a Hamiltonian differential system. The reason is that these equations can be expressed in terms of the Hamiltonian function H as [partial derivative][E.bar]/[partial derivative]t = [delta]H/[delta][H.bar] and [partial derivative][H.bar]/[partial derivative]H = -[[delta]H/[partial derivative][E.bar]]. In other words, the ideal Maxwell's equation have the following Hamiltonian functional H

H([H.bar],[E.bar]) = 1/2[[integral].sub.[OMEGA]]([1/[epsilon]][H.bar] x [nabla] x [H.bar] + [1/[mu]][E.bar] x [nabla] x [E.bar])d[OMEGA]. (14)

It is also worthy to point out here that Equations (2)-(3) conserve the following two invariants [11]

[W.sub.1](t) = [[integral].sub.[OMEGA]]([epsilon] [E.bar] x [E.bar] + [mu] [H.bar] x [H.bar]) d[OMEGA] (15)

[W.sub.2](t) = [[integral].sub.[OMEGA]]([epsilon][[absolute value of [partial derivative][E.bar]/[partial derivative]t].sup.2] + [mu][[absolute value of [partial derivative][H.bar]/[partial derivative]t].sup.2]) d[OMEGA]. (16)

Note that [W.sub.1](t) and [W.sub.2](t) are defined as the energy density I and the energy density II. When solving the canonical equations in (13), they should be integrated in the physical domain [OMEGA] using the symplectic type of integrators so as to retain the long-term accurate solution behavior.

To preserve symplectic structure and conserve total energy in the frequency independent Hamiltonian Maxwell's equations, the symplectic method of an implicit or an explicit type should be chosen. In addition to preserving the symplectic structure along the time direction, we are also aimed to develop a dispersively more accurate scheme in the spatial domain. Our strategy of getting a higher spatial accuracy is to make the numerical dispersion relation equation for the Maxwell's equations more closer to their exact dispersion relation equation. In the wavenumber space, it is therefore necessary to derive the numerical angular frequency explicitly in terms of the wavenumbers in x and y directions. Note that it is difficult or almost impossible to apply any symplectic type of the Runge-Kutta schemes, which involve several coupled solution steps, to derive the numerical dispersion relation equation. In this study the explicit symplectic Runge-Kutta scheme is employed so that it is possible to minimize the difference between the exact and numerical dispersion relation equations for the separable Hamiltonian system of Maxwell's equations.

It has been well known for some time that Maxwell's equations can be rewritten to the infinite dimensional Hamiltonian system which is endowed with the Hamiltonian H. We can therefore write the separable Hamiltonian function H (or the energy density) as the sum of T([E.bar]) and V([H.bar]), where T and V denote the kinetic energy and the potential energy, respectively. Since the Hamiltonian system of Maxwell's equations is separable, the explicit symplectic partitioned Runge-Kutta time-stepping scheme can be employed to integrate this differential system of equations through the tableaus for the respective Faraday's and Ampere's equations.

In this study, the second-order accurate explicit partitioned Runge-Kutta scheme suitable for the currently investigated separable Hamiltonian differential system is adopted to approximate the time derivative terms. The resulting semi-discretized equations are as follows

[[H.bar].sup.*] = [[H.bar].sup.n] - [dt/2[mu]][nabla] x [[E.bar].sup.n], (17)

[[E.bar].sup.n+1] = [[E.bar].sup.n] + [dt/[epsilon]][nabla] x [[H.bar].sup.*], (18)

[[H.bar].sup.n+1] = [[H.bar].sup.*] - [dt/2[mu]][nabla] x [[E.bar].sup.n+l]. (19)

5. DISPERSION RELATION EQUATION PRESERVING SPATIAL DISCRETIZATION METHOD

5.1. Discretization of One-dimensional Ampere's and Faraday's Equations

Since the first-order Ampere's and Faraday's equations can be transformed respectively to their equivalent second-order equations, which are [[partial derivative].sup.2][E.bar]/[partial derivative][t.sup.2] = [1/[epsilon][mu]][[nabla].sup.2][E.bar] and [[partial derivative].sup.2][H.bar]/[partial derivative][t.sup.2] = 1/[epsilon][mu]][[nabla].sup.2][H.bar], the spatial derivative terms for [E.bar] and [H.bar], which are [partial derivative][E.sub.z]/[partial derivative]x and [partial derivative][H.sub.y]/[partial derivative]x in the currently investigated TM-mode equations, should be approximated by a center-type numerical scheme. We first approximate the semi-discretized one-dimensional Faraday's equation by [[H.bar].sup.*] = [[H.bar].sub.n] - [dt/2[mu]][nabla] x [[E.bar].sup.n] using the physically correct center spatial scheme for the term [partial derivative][E.sub.z]/[partial derivative]x. At an interior point i in a grid system of uniform mesh size [DELTA]x, our approximation of the term [partial derivative][E.sub.z]/[partial derivative]x given below involves using six stencil points at i [+ or -] 1/2, i [+ or -] 3/2, i [+ or -] 5/2

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

The other first-order spatial derivative terms can be also expressed by the form shown in (20).

By substituting first the centered approximation equation for [[partial derivative][E.sub.z]/[partial derivative]x[|.sup.n] into Equation (17) and then substituting the resulting magnetic field solution [[H.bar].sup.*] into Equation (18), the discretized equation given below can be obtained for the one-dimensional Ampere's equation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (21)

Derivation of the above equation involves using [[H.bar].sup.n] = [[H.bar].sup.0] - [dt/2[mu]][nabla] x [[E.bar].sup.n] and [nabla] x [[H.bar].sup.0] = [[epsilon]/[DELTA]t]([[E.bar].sup.n] - [[E.bar].sup.n-1]), where [H.sup.0] is denotes the value of [H.bar] at t = (n - 1/2)[DELTA]t.

Following the same discretization procedures, the algebraic equation for the Faraday's equation can be derived as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (22)

As before, the above equation is derived using the equations given by [[H.bar].sup.n] = [[H.bar].sup.0] - [[DELTA]t/2[mu]][nabla] x [[E.bar].sup.n], [[E.bar].sup.n] = [[E.bar].sup.n-1] + [[DELTA]t/[epsilon]][nabla] x [[H.bar].sup.0], and [[H.bar].sup.0] = [[H.bar].sup.n-1] - [[DELTA]t/2[mu]][nabla] x [[E.bar].sup.n-1].

The free parameters [a.sub.1], [a.sub.2] and [a.sub.3] will be determined below to close the algebraic system of the derived discrete Maxwell's equations. To rigorously derive the algebraic equations for [a.sub.1], [a.sub.2] and [a.sub.3] for getting a good overall accuracy, our strategy is to reduce the amplitude as well as the phase errors generated from the discretized Ampere's and Faraday's equations. The modified equation analysis of second kind is performed first on Equation (21) or (22) by expanding the terms 0i[+ or -]5/2, [[phi].sub.i[+ or -]3/2] and [[phi].sub.i[+ or -]1/2], where [phi] = E or H, in Taylor series with respect to 0i. The two leading discretization errors shown in the resulting derived dispersion relation equation are then eliminated to get the following two equations

25[a.sup.2.sub.1] + 9[a.sup.2.sub.2] + [a.sup.2.sub.3] + 30[a.sub.1][a.sub.2] + 10[a.sub.1][a.sub.3] + 6[a.sub.2][a.sub.3] + 5[a.sub.1] + 3[a.sub.2] + [a.sub.3] = 0, (23)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (24)

To get the propagation characteristics of the TM-mode Maxwell's equations, we need, in particular, to reduce the dispersive error. The reason is that such an error can cause the phase velocity to be the functions of frequency and propagation angle. Any small error generated in the wave propagation velocity may be cumulated to reach an unacceptable level. The consequence is that the long-term simulation quality can be seriously deteriorated.

Dispersion relation equation is a good measure to examine how the angular frequency [omega] of the wave will vary with the wavenumber k. A higher dispersion accuracy can be obtained provided that the numerical angular frequency for the differential system of Ampere's and Faraday's equations relates well with the wavenumber. To derive the last algebraic equation, we are aimed to reduce the difference between the numerical dispersion relation equation and the exact dispersion relation equation. To achieve this goal, the equation in the space-time domain (x, t) is transformed to its corresponding angular frequencywavenumber space ([omega],kx). By substituting first the following planewave solution for [E.sub.z], which is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and the plane wave solution for [H.sub.y], which is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], into the differential Equation (21) or (22), the equation relating the angular frequency [omega] with the wavenumber [k.sub.x] can be derived as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (25)

For the purpose of minimizing the dispersion error, the equation relating the angular frequency [[omega].sub.exact] analytically with the wavenumber [k.sub.x] given below is employed for the one-dimensional Maxwell's equations

[[omega].sup.2.sub.exact] = [k.sup.2.sub.x][c.sup.2]. (26)

For getting a dispersive accuracy as high as possible, in this study the difference between the exact group velocity [partial derivative][[omega].sub.exact]/[partial derivative][k.sub.x] and the numerical group velocity [partial derivative][[omega].sub.num]/[partial derivative][k.sub.x] is minimized. Following the line of this thought, the error function defined as [[[partial derivative][[omega].sub.num]/[partial derivative][k.sub.x] - [partial derivative][[omega].sub.exact]/[partial derivative][k.sub.x]].sup.2] is minimized in a weak sense within the integral range given below

E = [[integral].sup.m[pi].sub.-m[pi]][[[partial derivative][[omega].sub.num]/[partial derivative][k.sub.x]] - [[partial derivative][[omega].sub.exact]/[partial derivative][k.sub.x].sup.2]W[gamma]d[gamma]. (27)

In the above, [gamma] = [k.sub.x] [DELTA]x denotes the scaled wavenumber and W([gamma]) is the weighting function. Inclusion of the weighting function in (27) is to make the integration of Equation (27) possible.

The error function E will be minimized by enforcing the limiting condition [partial derivative]E/[partial derivative][a.sub.3] = 0 to get the third algebraic equation. The equation derived from the above minimization means will be used together with the other two algebraic equations derived previously by the modified equation analysis of second kind. According to the simulation results, the best result is obtained at m = [pi]/2. The resulting three introduced coefficients [a.sub.i] (i = 1 ~ 3) in Equation (20) are [a.sub.1] = -0.00805, [a.sub.2] = 0.07443 and [a.sub.3] = -1.18302. Through the minimization study performed in the wavenumber domain and the modified equation analysis of second kind for[partial derivative][H.sub.x]/[partial derivative]x, the proposed center scheme with the best numerical dispersion relation equation is known to have the spatial accuracy order of fourth since [partial derivative][H.sub.x]/[partial derivative]x = [partial derivative][H.sub.x]/[partial derivative]x[|.sub.exact] - 0.015115 [h.sup.4][[[partial derivative].sup.5][H.sub.x]/[partial derivative][x.sup.5] + O([h.sup.6]) + ....

For completeness, the numerical angular frequency [[omega].sub.num] derived in this study is plotted against the wavenumber [kappa] in Figs. 1 and 2 at different values of Cr and m, respectively. Furthermore, we can see from Fig. 3 that our results are better than those predicted by the Yee's scheme.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

5.2. Discretization of Two-dimensional Ampere's and Faraday's Equations

The time derivative term in the two-dimensional equation is approximated first by the partitioned Runge-Kutta scheme to get the semi-discretized Faraday's equation. We then apply as before the physically correct center scheme to approximate the spatial terms [partial derivative][E.sub.z]/[partial derivative]x and [partial derivative][E.sub.z]/[partial derivative]y at an interior point (i,j) in a uniform grid system of the mesh size [DELTA]x = [DELTA]y. This center scheme development involves using the twelve nodal solutions at the stencil points (i [+ or -] 1/2, j), (i [+ or -] 3/2, j), (i [+ or -] 5/2, j), (i, j [+ or -] 1/2), (i, j [+ or -] 3/2) and (i, j [+ or -] 5/2) as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

Other first-order spatial derivative terms [partial derivative][H.sub.y]/[partial derivative]x and [partial derivative][H.sub.x]/[partial derivative]y are similarly expressed by the same center scheme as those shown in (28).

Substitution of the above centered equation for [partial derivative][E.sub.z]/[partial derivative]x[|.sup.n] into Equation (17) and then the magnetic field solution [[H.bar].sup.*] into Equation (18) enables us to get the following equation thanks to [[H.bar].sup.n] = [[H.bar].sup.0] - [[DELTA]t/2[mu]][nabla] x [[E.bar].sup.n] and [nabla] x [[H.bar].sup.0] = [[epsilon]/[DELTA]t]([[E.bar].sup.n] - [[E.bar].sup.n-1]).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (29)

The numerical dispersion relation equation given below for the two-dimensional Ampere's and Faraday's equations can be similarly derived by substituting the plane wave solutions [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] into the discrete Equation (29)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (30)

The same procedure as that employed in Subsection 3.1 is adopted to minimize the difference between the two-dimensional exact and numerical dispersion relation equations

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (31)

The above real-valued error function has the minimum value, which is zero, provided that [partial derivative]E/[partial derivative][a.sub.3] = 0. Application of this minimization means enables us to get one more algebraic equation to close the algebraic system in the determination of the introduced parameters, which are [a.sub.1] = 0.022277, [a.sub.2] = --0.17244 and [a.sub.3] = 0.77805, derived at Cr = 0.2 and m = [pi]/2. The numerical angular frequency [[omega].sub.num] is plotted against different flow angles in Fig. 4. For completeness, the group velocity is also plotted against the wavenumber [kappa] in Fig. 5 at different angles. For the sake of comparison, the numerical angular frequencies w are plotted against the wavenumber [kappa] in Fig. 6 using different schemes. Moreover, we compare the numerical dispersion errors predicted by different schemes in Fig. 7, from which our scheme is clearly seen to have a higher dispersive accuracy.

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

Given the definition of [k.sup.2] = [k.sup.2.sub.x] + [k.sup.2.sub.y], [k.sub.x] and [k.sub.y] can be expressed in terms of k and [theta] as [k.sub.x] = k cos [theta] and [k.sub.y] = k sin [theta]. For the sake of comparison of the numerical schemes and the discussion of the predicted results, the parameter [N.sub.[lambda]] = [lambda]/h, which represents the number of grid points per wavelength [lambda](= 2[pi]/k), and Cr = c[DELTA]t/h, which is known as the CFL number, are defined. The speed of light c is chosen as the reference speed and h = [DELTA]x = [DELTA]y as the uniform grid size. Given the above two definitions, the numerical phase velocity [V.sub.p], which is defined as the ratio of the numerical angular frequency and the wavenumber k, can be derived. We can express [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] by virtue of Equation (30) to get [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Define [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and, then, the ratio of the numerical phase velocity [V.sub.p] = [[omega].sub.R]/k and the exact phase velocity c can be derived as [12]

[V.sub.p]/c = [[omega].sub.R]/ck = [[N.sub.[lambda]]/2[pi]Cr][tan.sup.-1])(I/R). (32)

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

For the sake of comparison, the ratio of the derived numerical and exact phase velocities [V.sub.p]/c for the three schemes investigated at different values of [N.sub.[lambda]] are plotted in Fig. 8. For the case with fewer grid points per wavelength, our scheme is clearly seen to improve the phase velocity. Also, all the schemes are seen to have a better performance in the region near [theta] = 45[degrees]. As [N.sub.[lambda]] increases, the numerical phase velocity asymptotically approaches the exact phase velocity. In Fig. 9, the derived free parameters [a.sub.1], [a.sub.2] and [a.sub.3] are plotted against [theta] for the current scheme implemented at different values of the Courant number Cr.

6. VALIDATION AND ASSESSMENT OF THE DUAL PRESERVING MAXWELL'S EQUATION SOLVER

6.1. Analytical Validation of the Maxwell's Equations in Free Space

The explicit symplectic scheme with the optimized numerical dispersion relation equation developed in staggered grids will be validated by solving the TM-mode Maxwell's equations at [mu] = 1 and [epsilon] = 1 in -1 [less than or equal to] x [less than or equal to] 1 and -1 [less than or equal to] y [less than or equal to] 1. In this study, the problem amenable to analytic solution given below will be solved subject to the initial solenoidal solutions [E.sub.z](x,y,0) = cos(-2t), [H.sub.x](x,y,0) = sin [theta] x cos(-2t) and [H.sub.y](x, y, 0) = - cos [theta] x cos(-2t)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (33)

According to the computed errors tabulated in Table 2, the proposed dual-preserving scheme in free space for the TM-mode Maxwell's equations is analytically validated.

Besides performing the above validation study, the Hamiltonian defined in (14) and the energy densities I and II given in (15)-(16) are calculated. Note that the Hamiltonian is trivially equal to zero in any two-dimensional TM-mode Maxwell's equations. As a result, only the predicted and exact energy densities are plotted against time in Fig. 10. One can clearly see from this figure that the computed [W.sub.1] and [W.sub.2] values do not change with time using the explicit partitioned Runge-Kutta symplectic scheme. The value of [nabla] x H is predicted to be almost equal to zero. The Gauss's law is satisfied discretely using the explicit partitioned Runge-Kutta symplectic scheme.

[FIGURE 11 OMITTED]

6.2. Validation of the Maxwell's Equations in the Domain

Analysis is then performed in the square of area 4.694E-03 [m.sup.2] containing the dispersive media and the CPML layer with 10 cells. At the square centroid the incident Gaussian pulse with the bandwidth of 5-30 GHz is given below

Gaussian pulse = cos (2[pi]ft) x exp [-4[pi][(t - [t.sub.0]).sup.2]/[d.sup.2]. (34)

In the above, d = 200[DELTA]t, [t.sub.0] = 6/f and f = 17.5 GHz. The calculation carried out in a uniform mesh of [81 x 81] nodal points involves the grid size [DELTA]x = [DELTA]y = 8.565E - 04m and time step [DELTA]t = [DELTA]x/5c, where c is the speed of light. One can observe from Fig. 11 that the computed reflective errors at point P, which is at (70[DELTA]x, 70[DELTA]y) location, and Q, which is at (70[DELTA]x, 41[DELTA]y) in three investigated dispersive media are all negligibly small.

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

6.3. Assessment of the Proposed and Referenced Numerical Methods

The physical domain of current interest involves the left half of the vacuum and the right half contains the dispersive media, which include the three investigated Debye, Lorentz and Drude materials. Each half of the domain is divided into the 250 uniform intervals of the length [DELTA]x = 7.5 x [10.sup.-5] m. At two ends of the one-dimensional domain, perfect matching layers of length 10[DELTA]x are attached so as to absorb the possibly reflected wave.

A Gaussian pulse placed in the middle of the vacuum, or located at x = 125[DELTA]x, takes the form exp[-(t-[[t.sub.0]/d).sup.2]]. This pulse is defined by d = 40[DELTA]t, [t.sub.0] = 120[DELTA]t, where [DELTA]t = 0.125 x [10.sup.-12] second. The parameters of the three dispersive media on the right-half of the physical domain are chosen to be [[epsilon].sub.[infinity]] = 7, [[epsilon].sub.s] = 10, [tau] = 7.0 x [10.sup.-12]s for the Debye medium, [[epsilon].sub.[infinity]] = 1, [f.sub.p] = 28.7 GHz, [[omega].sub.p] = 2[pi][f.sub.p], [[gamma].sub.p] = 20 GHz for the Drude medium, and [[epsilon].sub.[infinity]] = 1.5, [[epsilon].sub.s] = 3, [[omega].sub.p] = 40[pi] GHz, [[delta].sub.p] = 4[pi] GHz for the Lorentz medium. The computed values of [E.sub.z] for all three dispersive media are plotted against time at the monitoring point 489[DELTA]x. One can clearly see that our results plotted in Fig. 12 agree well with the referenced solutions [13], thereby validating again the proposed numerical method.

With the success of validating the method for the one-dimensional simulation in dispersive medium, we will then apply the proposed method to simulate the problem involving the same Gaussian source as that in the 1D problem. In the physical domain of area 6.806E-05m2, the dispersive media of Debye, Lorentz and Drude types are surrounded by the perfect matching layer with the width of 10[DELTA]x to effectively absorb the outgoing wave. The simulated time-varying results of [E.sub.z] at the monitoring point, which is located at (100[DELTA]x, 100[DELTA]y), for the optical media of simple, Debye, Lorentz and Drude types are plotted in Fig. 13. As this figure shows, agreement between the present and the referenced solutions [13] is excellent.

7. NUMERICAL RESULTS IN DISPERSIVE MEDIA

7.1. Numerical Results in Debye Medium

We considered a one-dimensional problem consisting of 1000 cells. The first five hundred cells are used to model the vacuum (air) while the rest of five hundred cells are used for the Debye model (water) [10]. Each cell has a length of 37.5 [micro]m and the time step was 0.0625 ps. The parameters used in this example are [[epsilon].sub.s] = 81, [[epsilon].sub.[infinity]] = 1.8, and [t.sub.0] = 9. 4E - 12 [14]. The simulated propagation of this pulse through the air-water domain is shown in Fig. 14. The incident Gaussian pulse with the frequency spectrum up to 100 GHz is considered and is plotted in Fig. 15. The reflection occurring at the air-water interface (cell 500) is clearly visible. As the pulse propagates through the water, both of the attenuation and dispersion are clearly shown.

The wide bandwidth reflection coefficient was determined for the air-water interface by calculating the incident and reflected field strengths versus time at a position that is one cell ahead of the air-water interface. Two computations were performed. The first one involves the free-space parameters for the entire sample space. The result of this calculation shows the incident pulse at the interface versus time. A second calculation contains the air-water domain and the same position in front of the interface is considered. The result of the second calculation shows the total field in the front of the interface. The reflected field is obtained by subtracting the result of the incident field (no water) from the total field with the presence of water interface.

[FIGURE 14 OMITTED]

[FIGURE 15 OMITTED]

[FIGURE 16 OMITTED]

These electric fields versus time were transformed to the frequency domain by means of the Fast Fourier Transform (FFT). The reflection coefficient at each frequency was calculated by dividing the transforms of the reflected field and incident field results. Our results obtained after 4000 time steps (or at t = 250 ps) are compared well with the exact analytical frequency domain result in Fig. 16.

[FIGURE 17 OMITTED]

[FIGURE 18 OMITTED]

7.2. Numerical Results in Lorentz Medium

For the one-dimensional Lorentz problem [15], the parameters used for the current FDTD calculation involve the cell size [DELTA]x = 250 [micro]m, time increment [DELTA]t = 0.0625 ps (the Courant number limit), 1000 cells, and 2000 time steps. The dispersive half-space is characterized by a single pair of complex conjugate poles with the parameters [[epsilon].sub.s] = 3.0, [[epsilon].sub.[infinity]] = 1.5, [[omega].sub.p] = 2[pi] x 20 GHz and [[delta].sub.p] = 0.1[[omega].sub.p]. The conductivity a is zero. The propagation of the pulse through the vacuum-Lorentz space is shown in Fig. 17. The incident Gaussian pulse used in the present FDTD solution has a width of 40 cells and contains a significant amount of energy at the frequency up to approximately 100 GHz in Fig. 18.

The reflection coefficient is calculated by computing the FFT result for the time history of the reflected pulse at the interface and the FFT result for the time history of the incident pulse. Each frequency component of the first FFT result is divided by the corresponding component of the second FFT result to give the refection coefficient as the function of frequency. The simulated magnitudes of the reflection coefficient are compared well with the analytical data in Fig. 19.

7.3. Numerical Results in Drude Medium

We considered a problem consisting of 800 cells in the one-dimensional domain which includes the Drude medium (300 ~ 500 grids) and the other is vacuum. Each cell has the length of 75 [micro]m and the time increment was 0.125 ps and the number of time steps is 9600. The parameters used in this Drude model (plasma) are [[epsilon].sub.[infinity]] = 1.0, [f.sub.p] = 28.7 GHz, [[omega].sub.p] = 2[pi][f.sub.p], [[gamma].sub.p] = 20 GHz.

[FIGURE 19 OMITTED]

The propagation of the investigated pulse through the vacuum-plasma-vacuum space is simulated and the results are shown in Fig. 20. The incident wave was a derivative Gaussian pulse, shown in Fig. 21, with the frequency spectrum up to 100 GHz and the center frequency is 50 GHz. The reflection occurring at the vacuum-plasma interface (cell 299) is clearly visible. The transmission that occurs at the plasma-vacuum interface (cell 501) is also clearly visible. The attenuation and dispersion are also clearly exhibited as the pulse propagates through the plasma.

The wide bandwidth reflection coefficient was determined for the vacuum-plasma interface by calculating the incident and reflected field strengths versus time at a position with one cell ahead of the vacuum-plasma interface. The wide bandwidth transmission coefficient was determined for the plasma-vacuum interface by calculating the transmission field strength versus time at a position with one cell behind the plasma-vacuum interface.

Three computations were performed. The first is for the free-space parameter values. The result of this calculation provided the incident pulse at the interface (at cell number 299) versus time. A second calculation contains the vacuum-plasma domain. The result is shown at the same position in front of the interface. The result of the second calculation revealed the total field in front of the interface. The reflected field was obtained by subtracting the incident field result (no plasma) from the total field result with the presence of the plasma interface. The third calculation contains the plasma-vacuum domain and the result is plotted at the position (at cell number 501) in the rear of the plasma-vacuum interface. The result of the third calculation provided the transmission field.

These electric fields versus time were transformed to the frequency domain via the fast Fourier transform. The reflection coefficient at each frequency was calculated by dividing the transforms of the reflected field and the incident field. The transmission coefficient at each frequency was calculated by transforming the transmission field and the incident field. We compare our results after the 9600 time step with the exact frequency domain result for both of the reflection and transmission coefficients shown in Fig. 22.

[FIGURE 20 OMITTED]

[FIGURE 21 OMITTED]

[FIGURE 22 OMITTED]

8. CONCLUDING REMARKS

A theoretically rigorous solver is proposed to simulate the Maxwell's equations in frequency-independence and frequency-dependence media. Our aim is to generalize the Maxwell's equations solver in a sense that the symplectic property is present and the numerical dispersion relation equation is optimized in free space. All the simulation results obtained from the proposed method in Debye, Lorentz and Drude media agree excellently with the analytical as well as the benchmark solutions.

REFERENCES

[1.] Nicolaides, R. A. and D. Q. Wang, "Helicity and variational principles for Maxwell's equations," Int. J. Electron., Vol. 54, 861-864, 1983.

[2.] Cockburn, B., F. Li, and C. W. Chi, "Locally divergence-free discontinuous Galerkin methods for the Maxwell equations," Journal of Computational Physics, Vol. 194, No. 2, 588-610, 2004.

[3.] Sheu, W. H., Y. W. Hung, M. H. Tsai, P. H. Chiu, and J. H. Li, "On the development of a triple-preserving Maxwell's equations solver in non-staggered grids," Int. J. Numer. Meth. Fluids, Vol. 63, 1328-1346, 2010.

[4.] Sheu W. H., L. Y. Liang, and J. H. Li, "Development of an explicit symplectic scheme that optimizes the dispersion-relation equation of the Maxwell's equations," Communications in Computational Physics, Vol. 13, No. 4, 1107-1133, 2013.

[5.] Yee, K. S., "Numerical solution of initial boundary value problem involving Maxwell's equations in isotropic meida," IEEE Transactions on Antenna Propagation, Vol. 4, No. 3, 302-307, 1966.

[6.] Roden, J. A. and S. D. Gedney, "Convolutional PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media," Microwave Optical Tech. Lett., Vol. 27, 334-339, 2000.

[7.] Berenger, J. P., "A perfectly matched layer for the absorption of electromagnetic waves," Journal of Computational Physics, Vol. 114, No. 2, 185-200, 1994.

[8.] Taflove, A. and S. C. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method, 3rd Edition, Artech House, Norwood, MA, 2005.

[9.] Wei, B., S. Q. Zhang, F. Wang, and D. Ge, "A novel UPML FDTD absorbing boundary condition for dispersive media, waves in random and complex media," Journal of Mathematical Physics, Vol. 20, No. 3, 511-527, 2010.

[10.] Luebbers, R. J., F. P. Huusberger, K. S. Kunz, R. B. Standler, and M. Schneider, "A frequency-dependent finite-difference time-domain formulation for dispersive materials," IEEE Transactions on Electromagnetic Compatibility, Vol. 32, No. 3, 222-227, 1990.

[11.] Anderson, N. and A. M. Arthurs, "Helicity and variational principles for Maxwell's equations," Int. J. Electron, Vol. 54, 861-864, 1983.

[12.] Gao, L., B. Zhang, and D. Liang, "The splitting finite-difference time-domain methods for Maxwell's equations in two dimensions," J. Comput. Applied Math, Vol. 205, 207-230, 2007.

[13.] Wei, B., X. Y. Li, F. Wang, and D. Ge, "A finite difference time domain absorbing boundary condition for general frequency-dispersive media," Acta Physica Sinica, Vol. 58, No. 7, 6174-6178, 2009.

[14.] Cole, K. S. and R. H. Cole, "Dispersion and absorption in dielectrics," J. Chem. Phys, Vol. 9, 341, 1941.

[15.] Kelley, F. and R. J. Luebber, "Piecewise linear recursive convolution for dispersive media using FDTD," IEEE Transactions on Electromagnetic Compatibility, Vol. 44, No. 6, 1996.

T. W. H. Sheu (1, 2 3), *, R. Y. Chung (1), and J. H. Li (1)

(1) Department of Engineering Science and Ocean Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei, Taiwan, Republic of China

(2) Taida Institute of Mathematical Science (TIMS), National Taiwan University, Taiwan, Republic of China

(3) Center of Advanced Studies in Theoretical Science (CASTS), National Taiwan University, Taiwan, Republic of China

* Corresponding author: Tony W. H. Sheu (twhsheu@ntu.edu.tw).

Received 9 August 2012, Accepted 25 September 2012, Scheduled 15 October 2012

Table 1. Summary of the three investigated optical media. In this table, [[alpha].sub.p] = [[beta].sub.p], P = 1, [[beta].sub.p] = [([[omega].sup.2.sub.p] - [[delta].sup.2.sub.p]).sup.1/2] and [[gamma].sub.p] = [[G.sub.p]/[[beta].sub.p]]([[epsilon].sub.s] - [[epsilon].sub.[infinity]])[[omega].sup.2.sub.p]. Medium type Susceptibility x(omega) Debye medium x([omega]) = [[epsilon].sub.s] - [[epsilon].sub.[omega]]/ 1 + j[omega][tau] Lorentz medium x([omega]) = [[epsilon].sub.s] - [[epsilon].sub.[omega]][[omega].sup.2.sub.p]/ [[omega].sup.2.sub.p] + 2j[omega][[delta].sub.p] - [[omega].sup.2] Drude medium x([omega]) = [[omega].sup.2.sub.p]/ [[omega].sup.2] - j[omega][[gamma].sub.p] Medium type Susceptibility x(t) Debye medium x(t) = [[[epsilon].sub.s] - [[Epsilon].sub.[omega]]/[tau]] [e.sup.-t/[tau]]u(t) Lorentz medium [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Drude medium [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Medium type Relative electric permittivity [[epsilon].sub.r]([omega]) Debye medium [[epsilon].sub.r]([omega]) = [[epsilon].sub.[omega]] + [[[epsilon]. sub.s] - [[epsilon].sub.[omega]]/ 1 + j[omega][tau] Lorentz medium [[epsilon].sub.r]([omega]) = [[epsilon].sub.[omega]] + [([[epsilon].sub.s] - [[epsilon].sub.[omega]])[[omega].sup.2.sub.p]/ [[omega].sup.2.sub.p] + 2j[omega][[delta].sub.p] - [[omega].sup.2] Drude medium [[epsilon].sub.r]([omega]) = [[epsilon].sub. [infinity]] - [[omega].sup.2.sub.p]/[[omega]. sup.2] - j[omega][[gamma].sub.p] Table 2. Comparison of the computed [L.sub.2]-error norms and the required CPU times needed to get the results at t = 50 (s). The results are obtained at Cr = 0.2 and 0[degrees] [less than or equal to] 6 [less than or equal to] 90[degrees]. The weighting coefficients for [a.sub.1], [a.sub.2] and [a.sub.3] obtained at different angles are also tabulated for the sake of completeness. Degree [E.sub.z] CPU [a.sub.1] [L.sub.2]-error time(s) norm 0&90 3.6512E-05 22.64 -0.0073316 6&84 3.0065E-05 23.17 -0.0075212 9&81 6.1265E-05 23.50 -0.0074934 12&78 2.1540E-05 22.89 -0.0073804 22.5&67.5 4.9418E-05 25.08 -0.0068729 30&60 1.8921E-05 20.39 -0.0063184 36&54 2.1185E-05 20.78 -0.0058114 45 5.4198E-05 22.37 -0.0053241 Degree [a.sub.2] [a.sub.3] 0&90 0.074991 -1.1883 6&84 0.075939 -1.1902 9&81 0.075800 -1.1899 12&78 0.075236 _1.1888 22.5&67.5 0.072698 -1.1837 30&60 0.069925 _1.1782 36&54 0.067390 _1.1731 45 0.064954 -1.1682

Printer friendly Cite/link Email Feedback | |

Author: | Sheu, T.W.H.; Chung, R.Y.; Li, J.H. |
---|---|

Publication: | Progress In Electromagnetics Research |

Article Type: | Report |

Geographic Code: | 9CHIN |

Date: | Nov 1, 2012 |

Words: | 8154 |

Previous Article: | A wideband dual-polarized patch antenna with electric probe and magnetic loop feeds. |

Next Article: | Simple formula and its exact analytic solution of mutual impedance for nonplanar-skew dipoles. |

Topics: |