A Path Integral Representation Model to Extend the Analytical Capability of the Nonstandard Finite-difference Time-domain Method

Tadao Ohtani${}^{\mathbf{\text{1}}}$, Yasushi Kanai${}^{\mathbf{\text{2}}}$, and Nikolaos V. Kantartzis${}^{\mathbf{\text{3}}}$

¹Independent Researcher, Asahikawa, 070–0841, Japan
bytcg100@ybb.ne.jp

²Department of Engineering, Faculty of Engineering
Niigata Institute of Technology, Kashiwazaki 945–1195, Japan
kanai@iee.niit.ac.jp

³Department of Electrical and Computer Engineering
Aristotle University of Thessaloniki, Thessaloniki GR-54124, Greece
kant@ece.auth.gr

Submitted On: May 18, 2023; Accepted On: October 22, 2023

ABSTRACT

The nonstandard finite-difference time-domain (NS-FDTD) method is a powerful tool for solving Maxwell’s equations in their differential form on orthogonal grids. Nonetheless, to precisely treat arbitrarily shaped objects, very fine lattices should be employed, which often lead to unduly computational requirements. Evidently, such an issue hinders the applicability of the technique in realistic problems. For its alleviation, a new path integral (PI) representation model, equivalent to the NS-FDTD concept, is introduced. The proposed model uses a pair of basic and complementary path integrals for the H-nodes. To guarantee the same accuracy and stability as the NS-FDTD method, the two path integrals are combined via optimization parameters, derived from the corresponding NS-FDTD formulae. Since in the PI model, E-field computations on the complementary path are not necessary, the complexity is greatly reduced. Numerical results from various real-world problems prove that the proposed method improves notably the efficiency of the NS-FDTD scheme, even on coarse orthogonal meshes.

Index Terms: Computational electromagnetics, finite-difference time-domain methods, integral equations, radar cross section.

I. INTRODUCTION

Since its initial advent, the nonstandard finite-difference time-domain (NS-FDTD) method remains a very accurate scheme for electromagnetic field problems [1–15]. Indeed, its accuracy is 10${}^4$ times higher than that of the FDTD technique [1–4], in the case of coarse grids. Thus, the NS-FDTD method can be suitable for the electromagnetic design and analysis of electrically large structures with various dielectric materials, such as aircraft. Explicitly, the dimensions of a typical aircraft are approximately 500$\lambda$-1800$\lambda$, where $\lambda$ is the radar wavelength. Since the wavelength error of the FDTD method, along the axial direction, is 7$\lambda$-25$\lambda$ [16–19], it is apparent that such designs are very demanding. On the other hand, the error of the NS-FDTD algorithm is zero [1–4]. Despite this advantage, however, the method is established in differential form and applied to discrete points on orthogonal grids. Thus, its modeling accuracy degrades in the case of real objects with curved surfaces not aligned to the grid axes [16–19].

To overcome these shortcomings, the simplest way is the use of very fine grids, yet at the expense of large overheads. To this aim, a contour-path (CP) model based on the integral form of Maxwell’s equations has been presented [16], [17], [20], [21]. As the path integral (PI) scheme, stemming from the Stokes theorem, can handle arbitrary shapes, it can be a potential candidate for such problems, even with coarse lattices. In this context, 2D and 3D PI models have been developed for the NS-FDTD method [22], [23]. Nevertheless, these models are complex for practical applications as they require E-node computations for both the complementary and basic paths. Moreover, the numerical stability condition [16], [19] of the prior PI model has not been elaborately derived. The optimization parameters for the PI calculations are approximate values, derived from the numerical dispersion equation. Thus, the efficiency of this PI model could be questionable in some scenarios, e.g., long-term stability and compatibility with the NS-FDTD method in real-world configurations.

In this paper, a new PI form, equivalent to the NS- FDTD formulation in two (2D) and three dimensions (3D), is developed. To achieve a high accuracy with the same isotropy and stability level as the NS-FDTD method, both the basic and complementary path integrals for the H-node calculations on square grids are enhanced. Unlike previous PI models [22], [23], the proposed one does not need the E-field calculation on the complementary path. The two PIs are combined by new optimization parameters, analytically derived from an equivalency requirement with the NS-FDTD formulae. So, an advanced PI representation model is derived, which can be promptly employed with the NS-FDTD technique. Numerical validations reveal that the featured model drastically extends the applicability of the NS-FDTD algorithm for curved objects on orthogonal grids, attaining high levels of accuracy and convergence.

Figure 1: (a) Basic and (b) complementary (red arrows) integral path on a square grid for the new 2D PI model.

II. 2D PATH INTEGRAL FORM

A. Formulation

Let us consider a PI model on a square lattice that offers the best accuracy for the NS-FDTD method. In the FDTD scheme, a typical propagation error arises along the 0${}^o$ and 45${}^o$ directions [16]. For this reason, we adopt two distinct integral paths - the basic and the 45${}^o$-rotated complementary one - that can mutually cancel their errors, as shown in Fig. 1. For the H${}_z$component, the integral form of Maxwell’s equations is given by

$\mu {\displaystyle \frac{\partial }{\partial t}}{\displaystyle {\int }_S}\mathbf{H}\cdot d\mathbf{S}=$	$-{\displaystyle {\int }_C}\mathbf{E}\cdot d\mathbf{l}\Rightarrow \mu {\displaystyle \frac{\partial {\overline{H}}_z(x,y)}{\partial t}}{\mathrm{\Delta }}^2$	(1)
$=$	$-[E_y(x+\delta ,y)-E_x(x,y+\delta )$
		$-E_y(x-\delta ,y)+E_x(x,y-\delta )]\mathrm{\Delta },$

for the basic path with $\delta$ $\equiv$ $\mathrm{\Delta }$/2, and

		$\mu {\displaystyle \frac{\partial {\tilde{H}}_z(x,y)}{\partial t}}2{\mathrm{\Delta }}^2$	(2)
		$=-[{\tilde{E}}_y(x+\delta ,y+\delta )-{\tilde{E}}_x(x-\delta ,y+\delta )$
		$-{\tilde{E}}_y(x-\delta ,y-\delta )+{\tilde{E}}_x(x+\delta ,y-\delta )]\sqrt{2}\mathrm{\Delta },$

for the complementary path. To guarantee the desired isotropy, the $\partial$$H$${}_z$/$\mathrm{\partial }$t terms, obtained from the integral forms of (1a) and (1b), are combined as

$$\frac{\partial H_z(x,y)}{\partial t}={\beta }_0^{2\mathrm{D}}\frac{\partial {\overline{H}}_z(x,y)}{\partial t}+(1-{\beta }_0^{2\mathrm{D}})\frac{\partial {\tilde{H}}_z(x,y)}{\partial t},$$

(2)

where ${\beta }_0^{2\mathrm{D}}$ is a 2D optimization parameter, to be derived below, that can finely tune the resulting explicit expressions to the same accuracy as those of the NS-FDTD method. Next, the E${}_y$ component, on the basic path, in (1a), is calculated as

$$\begin{array}{cc}& \epsilon \frac{\partial }{\partial t}{\int }_S\mathbf{E}\cdot d\mathbf{S}={\int }_C\mathbf{H}\cdot d\mathbf{l}\hfill \\ & \Rightarrow \epsilon \frac{\partial E_z(x+\delta ,y)}{\partial t}=\frac{H_z(x,y)-H_z(x+\mathrm{\Delta },y)}{\mathrm{\Delta }}\hfill \end{array},$$

(3a)

		$E_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y)=E_y^{t-\mathrm{\Delta }t/2}(x+\delta ,y)$	(3b)
		$+{\displaystyle \frac{s_{\omega }(\mathrm{\Delta }t)}{\epsilon s_k(\mathrm{\Delta })}}\left[H_z^t(x,y)-H_z^t(x+\mathrm{\Delta },y)\right].$

Note that in (1)-(3), the cell width $\mathrm{\Delta }$ and temporal derivatives $\partial f/\partial t$ are substituted with $s_k(\mathrm{\Delta })=2\sin (k\mathrm{\Delta }/2)/k$ and $\left(f^{t+\mathrm{\Delta }t}-f^t\right)/s_{\omega }(\mathrm{\Delta }t)$, respectively, where $k$ is the wave number, $\omega$ the angular frequency, $\mathrm{\Delta }t$ the time increment, and $s_{\omega }(\mathrm{\Delta }t)=2\sin (\omega \mathrm{\Delta }t/2)/\omega$ [1–3]. A similar analysis holds for the $E_x$ case, as well. Moreover, the ${\tilde{E}}_y$ component on the complementary path in (1b), is derived by its existing value on the basic path as

$$\begin{array}{c}\hfill {\tilde{E}}_y(x+\delta ,y+\delta )=\{[E_x(x,y+\delta )+E_x(x+\mathrm{\Delta },y+\delta )]\frac{{\mathbf{e}}_x}{2}\hfill \\ \hfill +[E_y(x+\delta ,y)+E_y(x+\delta ,y+\mathrm{\Delta })]\frac{{\mathbf{e}}_y}{2}\}\cdot \frac{-{\mathbf{e}}_x+{\mathbf{e}}_y}{\sqrt{2}},\hfill \end{array}$$

(4)

while an analogous set for the E${}_x$ and ${\tilde{E}}_x$ components can be equivalently acquired. Also, in (4), e${}_x$${}_{\mathrm{,}}$${}_y$ are the unit vectors on the basic path and $(-{\mathbf{\text{e}}}_x+{\mathbf{\text{e}}}_y)/\sqrt{2}$ the unit vector on the complementary path. Therefore, for the H${}_z$ component, in (2), we obtain

$H_z^{t+\mathrm{\Delta }t}(x,y)=H_z^t(x,y)$
$-{\beta }_0^{2\mathrm{D}}{\displaystyle \frac{s_{\omega }(\mathrm{\Delta }t)}{\mu s_k(\mathrm{\Delta })}}[E_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y)-E_x^{t+\mathrm{\Delta }t/2}(x,y+\delta )$
$-E_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y)+E_x^{t+\mathrm{\Delta }t/2}(x,y-\delta )]$
$-(1-{\beta }_0^{2\mathrm{D}}){\displaystyle \frac{s_{\omega }(\mathrm{\Delta }t)}{\mu \sqrt{2}{s}_k(\mathrm{\Delta })}}[{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y+\delta )$
$-{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x-\delta ,y+\delta )-{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y-\delta )$
$+{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x+\delta ,y-\delta )],$	(5)

which is the formula for the H${}_z$ node of our PI model. It is stressed that the use of the interpolated ${\tilde{E}}_y$ value, from (4), in the PI scheme is different from the 2D PI model of [22], due to the complementary path integral.

B. Derivation of the ${\beta }_0^{2\mathrm{D}}$ parameter

To obtain the optimal ${\beta }_0^{2\mathrm{D}}$ parameter, we compare the PI form of (2) and the NS-FDTD formula with the operators [1–3]

$$d_x^{(0)}={\alpha }_0{d}_x^{(1)}+\left(1-{\alpha }_0\right)d_x^{(2)},$$

(6a)

$$d_x^{(1)}f(x,y)=f(x+\delta ,y)-f(x-\delta ,y),$$

(6b)

$d_x^{(2)}f(x,y)=[f(x+\delta ,y-\mathrm{\Delta })-f(x-\delta ,y-\mathrm{\Delta })$
$\mathrm{\hspace{1em}}+f(x+\delta ,y+\mathrm{\Delta })-f(x-\delta ,y+\mathrm{\Delta })]/2,$	(6c)

with analogous forms holding for the $d_y^{(*)}$ ones. Note that, herein, we focus on the $\partial$${}_x$E${}_y$ terms to acquire a simple argument. Hence, and after some calculus in (2), one may reach to

		$-{\mu {\displaystyle \frac{\partial H_z(x,y)}{\partial t}}|}_{{\partial }_x{E}_y}$	(7)
		$={\displaystyle \frac{1+{\beta }_0^{2\mathrm{D}}}{2}}\left[E_y(x+\delta ,y)-E_y(x-\delta ,y)\right]$
		$+{\displaystyle \frac{1-{\beta }_0^{2\mathrm{D}}}{4}}[E_y(x+\delta ,y+\mathrm{\Delta })-E_y(x-\delta ,y+\mathrm{\Delta })$
		$+E_y(x+\delta ,y+\mathrm{\Delta })-E_y(x-\delta ,y+\mathrm{\Delta })\left]\right|{}_{\mathrm{PI}}$
		$\equiv {\alpha }_0{d}_x^{(1)}{E}_y(x,y)+{\left(1-{\alpha }_0\right)d_x^{(2)}{E}_y(x,y)|}_{\mathrm{NS}-\mathrm{FDTD}},$

Comparing ${\beta }_0^{2\mathrm{D}}$ and $\alpha$${}_0$ from [1–3], we find that

$$\frac{1+{\beta }_0^{2\mathrm{D}}}{2}\equiv {\alpha }_0=\frac{1+{\gamma }_0}{2}\Rightarrow {\beta }_0^{2\mathrm{D}}={\gamma }_0=\frac{2}{3}-\frac{{(k\mathrm{\Delta })}^2}{90}.$$

(8)

Therefore, it is proven that (2.1) can be equated with the 2D NS-FDTD formula through ${\beta }_0^{2\mathrm{D}}$ given in (8). This implies that the wave propagation characteristics, accuracy, and stability of the new PI model are the same as those of the NS-FDTD method [1–3], [7], [13].

III. 3D PATH INTEGRAL FORM

A. Formulation

For the improvement of the 2D PI model, shown in Fig. 1, we must eliminate the propagation error due to angular dependence. To this objective, six square integral paths are devised; i.e., the basic and complementary paths are mutually rotated 45${}^o$ (around the z-axis) at the z $\pm$ $\mathrm{\Delta }$ planes, as in Fig. 2. Using basic and complementary paths in Fig. 2, our PI forms for the H${}_z$ component are

		$\mu {\displaystyle \frac{\partial }{\partial t}}{\displaystyle {\int }_S}\mathbf{H}\cdot d\mathbf{S}=-{\displaystyle {\int }_C}\mathbf{E}\cdot d\mathbf{l}\Rightarrow \mu {\displaystyle \frac{\partial {\overline{H}}_z(x,y,z)}{\partial t}}{\mathrm{\Delta }}^2$	(9a)
		$=-[E_y(x+\delta ,y,z)-E_x(x,y+\delta ,z)$
		$-E_y(x-\delta ,y,z)+E_x(x,y-\delta ,z)]\mathrm{\Delta },$

for the basic path, and

Figure 2: Basic and complementary integral paths for the calculation of the H${}_z$ component on a 3D square grid. Six paths are used at the z and z $\pm$ $\mathrm{\Delta }$ planes, while in the two figures on the right, H${}_z$(z) is, actually, H${}_z$(x, y, z) and $H$${}_z$(z $\pm$ $\mathrm{\Delta }$) is $H$${}_z$(x, y, z $\pm$ $\mathrm{\Delta }$).

$\mu$	${\displaystyle \frac{\partial {\tilde{H}}_z(x,y,z)}{\partial t}}2{\mathrm{\Delta }}^2$	(9b)
$=$	$-[{\tilde{E}}_y(x+\delta ,y+\delta ,z)-{\tilde{E}}_x(x-\delta ,y+\delta ,z)$
		$-{\tilde{E}}_y(x-\delta ,y-\delta ,z)+{\tilde{E}}_x(x+\delta ,y-\delta ,z)]\sqrt{2}\mathrm{\Delta },$

for the complementary paths at the z and z $\pm$ $\mathrm{\Delta }$ planes. To attain the desired isotropy, we combine the integral forms of (9a) and (9b) as

$$\begin{array}{c}\hfill \frac{\partial H_z(x,y,z)}{\partial t}={\beta }_1\left[{\beta }_0\frac{\partial {\overline{H}}_z(x,y,z)}{\partial t}+\left(1-{\beta }_0\right)\frac{\partial {\tilde{H}}_z(x,y,z)}{\partial t}\right]\\ \hfill +\frac{1-{\beta }_1}{2}[{\beta }_2\frac{\partial {\overline{H}}_z(x,y,z+\mathrm{\Delta })}{\partial t}+(1-{\beta }_2)\frac{\partial {\tilde{H}}_z(x,y,z+\mathrm{\Delta })}{\partial t}\\ \hfill +{\beta }_2\frac{\partial {\overline{H}}_z(x,y,z-\mathrm{\Delta })}{\partial t}+(1-{\beta }_2)\frac{\partial {\tilde{H}}_z(x,y,z-\mathrm{\Delta })}{\partial t}]\end{array},$$

(10)

where $\beta$${}_{0,1}$ are the optimization parameters for an isotropic wave propagation. Basically, (10) denotes our 3D PI model. Then, the E${}_y$ component in (9a) becomes

		$\epsilon {\displaystyle \frac{\partial }{\partial t}}{\displaystyle {\int }_S}\mathbf{E}\cdot d\mathbf{S}={\displaystyle {\int }_C}\mathbf{H}\cdot d\mathbf{l}\Rightarrow \epsilon {\displaystyle \frac{\partial E_y(x+\delta ,y,z)}{\partial t}}{\mathrm{\Delta }}^2$	(11a)
		$=[H_z(x,y,z)-H_z(x+\mathrm{\Delta },y,z)$
		$+H_x(x+\delta ,y,z+\delta )-H_x(x+\delta ,y,z-\delta )]\mathrm{\Delta },$

		$E_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y,z)=E_y^{t-\mathrm{\Delta }t/2}(x+\delta ,y,z)$	(11b)
		$+{\displaystyle \frac{s_{\omega }(\mathrm{\Delta }t)}{\epsilon s_k(\mathrm{\Delta })}}[H_z^t(x,y,z)-H_z^t(x+\mathrm{\Delta },y,z)$
		$+H_x^t(x+\delta ,y,z+\delta )-H_x^t(x+\delta ,y,z-\delta )].$

		$H_z^{t+\mathrm{\Delta }t}(x,y,z)=H_z^t(x,y,z)$	(13)
		$-{\displaystyle \frac{{\beta }_1{s}_{\omega }(\mathrm{\Delta }t)}{\mu s_k(\mathrm{\Delta })}}\{{\beta }_0[E_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y,z)-E_x^{t+\mathrm{\Delta }t/2}(x,y+\delta ,z)-E_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y,z)+E_x^{t+\mathrm{\Delta }t/2}(x,y-\delta ,z)]$
		$+{\displaystyle \frac{1-{\beta }_0}{\sqrt{2}}}[{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y+\delta ,z)-{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x-\delta ,y+\delta ,z)-{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y-\delta ,z)+{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x+\delta ,y-\delta ,z)]\}$
		$-{\displaystyle \frac{\left(1-{\beta }_1\right)s_{\omega }(\mathrm{\Delta }t)}{2\mu s_k(\mathrm{\Delta })}}\{{\beta }_2[E_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y,z+\mathrm{\Delta })-E_x^{t+\mathrm{\Delta }t/2}(x,y+\delta ,z+\mathrm{\Delta })-E_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y,z+\mathrm{\Delta })+E_x^{t+\mathrm{\Delta }t/2}(x,y-\delta ,z+\mathrm{\Delta })]$
		$+{\displaystyle \frac{1-{\beta }_2}{\sqrt{2}}}[{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y+\delta ,z+\mathrm{\Delta })-{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x-\delta ,y+\delta ,z+\mathrm{\Delta })-{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y-\delta ,z+\mathrm{\Delta })+{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x-\delta ,y-\delta ,z+\mathrm{\Delta })]\}$
		$-{\displaystyle \frac{\left(1-{\beta }_1\right)s_{\omega }(\mathrm{\Delta }t)}{2\mu s_k(\mathrm{\Delta })}}\{{\beta }_2[E_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y,z-\mathrm{\Delta })-E_x^{t+\mathrm{\Delta }t/2}(x,y+\delta ,z-\mathrm{\Delta })-E_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y,z-\mathrm{\Delta })+E_x^{t+\mathrm{\Delta }t/2}(x,y-\delta ,z-\mathrm{\Delta })]$
		$+{\displaystyle \frac{1-{\beta }_2}{\sqrt{2}}}[{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x+\delta ,y+\delta ,z-\mathrm{\Delta })-{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x-\delta ,y+\delta ,z-\mathrm{\Delta })-{\tilde{E}}_y^{t+\mathrm{\Delta }t/2}(x-\delta ,y-\delta ,z-\mathrm{\Delta })+{\tilde{E}}_x^{t+\mathrm{\Delta }t/2}(x-\delta ,y-\delta ,z-\mathrm{\Delta })]\}.$

		$-{\mu \mathrm{\Delta }{\displaystyle \frac{\partial H_z(x,y,z)}{\partial t}}|}_{\partial x}={\displaystyle \frac{{\beta }_1\left(1+{\beta }_0\right)}{2}}\left[E_y(x+\delta ,y,z)-E_y(x-\delta ,y,z)\right]$	(15)
		$+{\displaystyle \frac{\left(1-{\beta }_1\right)\left(1-{\beta }_2\right)}{8}}[E_y(x+\delta ,y+\mathrm{\Delta },z+\mathrm{\Delta })-E_y(x-\delta ,y+\mathrm{\Delta },z+\mathrm{\Delta })+E_y(x+\delta ,y-\mathrm{\Delta },z+\mathrm{\Delta })-E_y(x-\delta ,y-\mathrm{\Delta },z+\mathrm{\Delta })$
		$+E_y(x+\delta ,y+\mathrm{\Delta },z-\mathrm{\Delta })-E_y(x-\delta ,y+\mathrm{\Delta },z-\mathrm{\Delta })+E_y(x+\delta ,y-\mathrm{\Delta },z-\mathrm{\Delta })-E_y(x-\delta ,y-\mathrm{\Delta },z-\mathrm{\Delta })]$
		$+{\displaystyle \frac{{\beta }_1\left(1-{\beta }_0\right)}{4}}\left[E_y(x+\delta ,y+\mathrm{\Delta },z)-E_y(x-\delta ,y+\mathrm{\Delta },z)+E_y(x+\delta ,y-\mathrm{\Delta },z)-E_y(x-\delta ,y-\mathrm{\Delta },z)\right]$
		$+{\displaystyle \frac{\left(1-{\beta }_1\right)\left(1+{\beta }_2\right)}{4}}\left[E_y(x+\delta ,y,z+\mathrm{\Delta })-E_y(x-\delta ,y,z+\mathrm{\Delta })+E_y(x+\delta ,y,z-\mathrm{\Delta })-E_y(x-\delta ,y,z-\mathrm{\Delta })\right].$

As with the 2D PI formula, $\mathrm{\Delta }$ is replaced with to $s_k(\mathrm{\Delta })$ and $\partial f/\partial t$ with $\left(f^{t+\mathrm{\Delta }t}-f^t\right)/s_{\omega }(\mathrm{\Delta }t)$, according to the NS-FDTD concept. A similar procedure can also be performed for the E${}_x$ component. In this context, the interpolated ${\tilde{E}}_{x,y}$ quantities in (9b) are derived from

		${\tilde{E}}_x(x-\delta ,y+\delta ,z)$	(12a)
		$=\{[E_x(x-\mathrm{\Delta },y+\delta ,y,z)+E_x(x,y+\delta ,z)]{\mathbf{e}}_x$
		$+[E_y(x-\delta ,y,z)+E_y(x-\delta ,y+\mathrm{\Delta },z)]{\mathbf{e}}_y\}\cdot {\displaystyle \frac{{\mathbf{e}}_x+{\mathbf{e}}_y}{2\sqrt{2}}},$

${\tilde{E}}_y$	$(x+\delta ,y+\delta ,z)$	(12b)
		$=\{[E_x(x,y+\delta ,y,z)+E_x(x+\mathrm{\Delta },y+\delta ,z)]{\mathbf{e}}_x$
		$+[E_y(x+\delta ,y,z)+E_y(x+\delta ,y+\mathrm{\Delta },z)]{\mathbf{e}}_y\}\cdot {\displaystyle \frac{-{\mathbf{e}}_x+{\mathbf{e}}_y}{2\sqrt{2}}}.$

Therefore, the final path integral expression for the H${}_z$ component is given by (13) (see top of the next page). Note that the process of (9)-(13) can be analogously applied to the extraction of the H${}_{x\mathrm{,}y}$ components.

B. Derivation of the ${\beta }_0$, ${\beta }_1$, and ${\beta }_2$ parameters

Let us decompose the temporal derivative on the left-hand side of (10) into a $\partial$x and a $\partial$y term, for a simple argument, in order to derive the $\beta$${}_{0,1,2}$ parameters; i.e.,

$$\frac{\partial H_z(x,y,z)}{\partial t}={\frac{\partial H_z(x,y,z)}{\partial t}|}_{\partial x}+{\frac{\partial H_z(x,y,z)}{\partial t}|}_{\partial y}.$$

(14)

Then, the $\partial$x part of (14) is derived via (15) (see top of this page). Comparing (15) with the 3D NS-FDTD operators

$$d_x^{(0)}={\alpha }_1{d}_x^{(1)}+{\alpha }_2{d}_x^{(2)}+{\alpha }_3{d}_x^{(3)},$$

(16)

and by means of Appendix $\mathrm{A}$ and [1–3], [13], we obtain the following coefficient relations

$${\beta }_1\left(1+{\beta }_0\right)=2{\alpha }_1,$$

(17a)

$$\left(1-{\beta }_1\right)\left(1-{\beta }_2\right)=2{\alpha }_2,$$

(17b)

$${\beta }_1\left(1-{\beta }_0\right)=\left(1-{\beta }_1\right)\left(1+{\beta }_2\right)={\alpha }_3.$$

(17c)

If the system of equations (17a)-(17c) is solved by means of the $\alpha$${}_1$ + $\alpha$${}_2$ + $\alpha$${}_3$ =1 constraint, the new optimization parameters $\beta$${}_{0,1,2}$ can be extracted as

$$(17a)\&(17c):{\beta }_1={\alpha }_1+0.5{\alpha }_3,$$

(18a)

$$(18a)\&(17b):{\beta }_2=(0.5{\alpha }_3-{\alpha }_2)/(0.5{\alpha }_3+{\alpha }_2),$$

(18b)

$$(18a)\&(17c):{\beta }_0=({\alpha }_1-0.5{\alpha }_3)/({\alpha }_1+0.5{\alpha }_3).$$

(18c)

In this manner, the featured PI formula, given in (13), becomes fully equivalent to the corresponding 3D NS-FDTD expression, through the optimization parameters of (18a)-(18c), on cubic grids. Consequently, the wave propagation characteristics along with the numerical accuracy and stability of the developed PI model, in a homogeneous space discretized into cubic cells, are the same as those of the 3D NS-FDTD method.

IV. NUMERICAL RESULTS

In this section, we modify and extend our PI model to successfully treat curved structures. Then, it is verified whether it can be utilized as the CP technique of the FDTD method. Also, by combining this PI model with the NS-FDTD concept, we prove the superior performance of the former in realistic applications. So it is shown that the novel PI scheme significantly enhances the efficiency of the NS-FDTD algorithm for arbitrary object shapes, without opting for unduly fine lattices.

A. Analysis of 2D applications

The featured PI scheme is, firstly, applied to the radar cross section (RCS) analysis of a perfect electric conductor (PEC) circular cylinder with a radius of $\lambda$/2. The modification of the PI model near the cylinder is shown in Fig. 3. In particular, the integral area for the computation of the H${}_z$ component is enclosed by two paths, i.e., the S${}_B$and the S${}_C$, except for the PEC cylinder region. Note that the treatment of the path integrals is based on the typical FDTD principles [16]. So, the PI form for the basic path in Fig. 3 (a) can be written as

$\mu {\displaystyle \frac{\partial }{\partial t}}{\displaystyle {\int }_S}\mathbf{H}\cdot d\mathbf{S}$	$=-{\displaystyle {\int }_C}\mathbf{E}\cdot d\mathbf{l}\Rightarrow \mu {\displaystyle \frac{\partial {\overline{H}}_z(x,y)}{\partial t}}{S}_B^{\prime }{\mathrm{\Delta }}^2$	(19)
		$=E_x(x,y+\delta )\mathrm{\Delta }-E_y(x+\delta ,y)\mathrm{\Delta }$
		$+E_y(x-\delta ,y)l_y\mathrm{\Delta }-E_x(x,y-\delta )l_x\mathrm{\Delta },$

where $S_B^{\prime }=S_B/{\mathrm{\Delta }}^2$ and $l_{x,y}=l_{\mathit{\text{path length}}}/\mathrm{\Delta }$. Similarly, we can derive the PI expression for the complementary path by means of $S_C$, as depicted in Fig. 3 (b). It is mentioned that during the evaluation of (19), $\mathrm{\Delta }$ must be replaced with $s_k(\mathrm{\Delta })=2\sin (k\mathrm{\Delta }/2)/k$. Furthermore, due to the lack of $H_z$ nodes in the PEC area, there are some $E_{xy}$ and ${\tilde{E}}_{x,y}$ components that can not be computed via (3b) and (4). In such a case, these quantities are derived through the extrapolation/projection of already calculated E-field values on the nearest-neighbor nodes. The remaining E- and H-field components are acquired in terms of (3)-(5).

Figure 3: Modified (a) basic and (b) complementary path for the proposed PI model. Dashed lines indicate the position of the pair paths, and PEC areas are shown in light blue. The integral path and the area for the calculation of the H${}_z$ component are shown in the right-hand side figures (except the PEC region).

Figure 4: RCS variation versus the observation angle $\theta$ for diverse grid sizes $\mathrm{\Delta }$. The PEC cylinder is modeled via the staircase FDTD method, and $\mathrm{\Delta }$t = T/15, T/85, T/145, T/170 are the temporal increments that correspond to $\mathrm{\Delta }$ = $\lambda$/10, $\lambda$/60, $\lambda$/100, $\lambda$/120, respectively.

Figure 5: RCS variation for the 2D PEC cylinder, computed via the proposed PI model and the staircase NS-FDTD method only (without PI treatment). The incidence angle of the TE plane wave is $\theta$${}^{inc}$ = 0${}^o$ and $\mathrm{\Delta }$ = $\lambda$/10, with $\lambda$ = 1 m. The reference solution is obtained in terms of the staircase FDTD technique with $\mathrm{\Delta }$ = $\lambda$/120.

Regarding the reference solution of our problem, we examine several FDTD simulations of the PEC cylinder, as illustrated in Fig. 4. Since convergence is attained for $\mathrm{\Delta }=\lambda /100$ and $\lambda$/120, we select the $\mathrm{\Delta }$ = $\lambda$/120 grid resolution with $\mathrm{\Delta }$t = T/170, for T = 2$\pi$/$\omega$, as our reference. The structure is illuminated by a TE incident plane, with an incidence angle of $\theta$${}^{inc}$ = 0${}^o$, and all scattered waves are obtained by the TF/SF formulation [16–18] applied in the NS-FDTD region. Moreover, our domain is discretized into cells of $\mathrm{\Delta }$ = $\lambda$/10 ($\lambda$ = 1 m) for the basic path and $\mathrm{\Delta }$ = $\lambda$/7 for the complementary path, while stability is guaranteed with a $\mathrm{\Delta }$t = T/15. The PI area has a size of 14$\mathrm{\Delta }$$\times$14$\mathrm{\Delta }$ (including the cylinder surface) and is surrounded by a 500$\mathrm{\Delta }$$\times$500$\mathrm{\Delta }$ NS-FDTD as well as a 20$\mathrm{\Delta }$-thick perfectly matched layer (PML) [16–18].

In this context, Fig. 5 compares the RCS of the PEC cylinder computed via the proposed PI model and the NS-FDTD method solely (i.e., without PI treatment). As observed, the proposed technique is in very good agreement with the reference solution, despite the fairly coarse resolution of $\mathrm{\Delta }$ = $\lambda$/10 and $\mathrm{\Delta }$ = $\lambda$/7 for the basic and complementary paths, respectively. Concerning the small ripple in the PI results, we presume that it is attributed to the $\mathrm{\Delta }$ = $\lambda$/7 size of the complementary path cells. Actually, this ripple disappears in the case of a $\mathrm{\Delta }$ = $\lambda$/12 mesh resolution, as verified in Figs. 7 and 9 (a) (see below). So, a viable choice for retaining the balance between an acceptable accuracy level and a limited computational burden can be $\mathrm{\Delta }$ $\le$ $\lambda$/12. In addition, Fig. 5 indicates that the NS-FDTD method, without the PI model, fails to deliver acceptable accuracy. In summary, our PI model achieved a considerable (120/10)${}^2$ = 144 times reduction of the memory overhead with regard to the number of cells, and a (120/10)${}^2$$\times$ (170/15) = 1632 times speed improvement, with regard to the number of time steps, compared to those of the reference FDTD ($\mathrm{\Delta }$ = $\lambda$/120, $\mathrm{\Delta }$t = T/170) solution. Thus, the new PI concept can greatly enhance the efficiency of the NS-FDTD formulation in the case of arbitrarily-curvedobjects.

Figure 6: Cross-section of a PEC wing-like model. All dimensions are in $\mathrm{\Delta }$ units, with $\mathrm{\Delta }$ = $\lambda$/12 and $\lambda$ = 1 m.

Figure 7: RCS variation for the 2D wing-like model illuminated by an incident plane wave with (a) $\theta$${}^{inc}$ = 0${}^o$, (b) $\theta$${}^{inc}$ = 18${}^o$, and (c) $\theta$${}^{inc}$ = 25${}^o$, computed through the staircase NS-FDTD method only (without PI treatment, $\mathrm{\Delta }$ = $\lambda$/12, $\lambda$ = 1 m), the proposed PI model ($\mathrm{\Delta }$ = $\lambda$/12), and the reference staircase FDTD solution with $\mathrm{\Delta }$ = $\lambda$/96.

Next, we consider the RCS analysis of the PEC wing-like model, described in Fig. 6, whose dimensions are measured in $\mathrm{\Delta }$ units, with $\mathrm{\Delta }$ = $\lambda$/12 and $\lambda$ = 1 m. For its illumination, we use an incident, $H_z^{inc}$, plane wave. Moreover, the space around the structure is modeled via the NS-FDTD method and the modified PI scheme is applied to the cells near the wing-like surface. The RCS variation for different incidence angles, $\theta$${}^{inc}$, is presented in Fig. 7. The simulations are performed through the novel PI model, the staircase NS-FDTD technique only (i.e., without any PI treatment), and the reference FDTD solution with $\mathrm{\Delta }$ = $\lambda$/96. From the results, it is deduced that, in all cases, the proposed PI model accomplishes the best accuracy and coincidence with the reference solution, unlike the typical staircaseimplementations.

The prior 2D applications have also led to some interesting observations. The first one refers to the consistency of the geometrical accuracy of the integral path length with its area, given in Fig. 3. Hence, no discrepancies can arise during the combined basic and complementary path calculations. The second remark is related to the precision in the calculation of the E and $\tilde{E}$ quantities on the modified paths near the surface of the scatterer. Such computations are conducted via an extrapolation/projection process, which uses the inner product of the unit e${}_x$${}_{\mathrm{,}}$${}_y$vectors with the nearest-neighbor E field components, as occurs in (12). This procedure is important for our PI evaluations and the geometrical accuracy of the modified PI cells,as well.

Figure 8: Staircase model of a 12$\lambda$-long PEC cylinder with a circular cross-section of radius $\lambda$, for $\mathrm{\Delta }$ = $\lambda$/12 and $\lambda$ = 1 m. The modified PI cells are applied only to the surface of the structure.

B. Analysis of 3D applications

To prove the efficacy of the developed PI scheme in 3D problems, we focus on the RCS study of a finite-length PEC cylinder and the investigation of a thin-wire antenna impedance using a thin wire model [16], [25], [26]. The radius and the length of the former structure are much larger than $\mathrm{\Delta }$, while the latter has a fine form with a curved surface that is smaller than $\mathrm{\Delta }$.

1) RCS study of a finite-length PEC cylinder

The staircase model (cubic cells with $\mathrm{\Delta }$ = $\lambda$/12 and $\lambda$ = 1 m) of the 12$\lambda$-long PEC cylinder, whose circular cross-section has a radius of $\lambda$, is illustrated in Fig. 8. The incident plane wave is $E_{\varphi }^{inc}$, with its angles ($\theta$, $\varphi$) expressed in polar coordinates. In our simulations, the PI scheme of Fig. 3 replaces the staircase approximation only at the smooth surface of the cylinder. Figure 9 gives the RCS results for the E${}_{\varphi }$ component, evaluated via several techniques. Note that, apart from the schemes already employed in the previous examples, we herein compare our results with the PI formulation of [23], where the interpolated $\tilde{E}$ values are obtained from the PI scheme instead of (12). Again, the outcomes of the novel PI model are pretty close to the reference solution and in satisfactory agreement with the PI model of [23]. The presence of some small discrepancies between the two PI models, at the low RCS levels, is mainly due to the different treatment of the E-field components on the modified paths. Indeed, the proposed PI model retrieves the required $\tilde{E}$ terms by means of (12), whereas that of [23] derives them directly from Ampère’s law. Nonetheless, another possible reason could be the different handling of those $\tilde{E}$ terms that cannot be promptly computed owing to the presence of the PEC object. In any case, our scheme reduces the overall burden, regarding the number of cells, to (72/12)${}^3$= 216 times, compared to the staircase FDTD method. Hence, the prior facts substantiate the competence of the model to successfully manipulate 3D curved structures.

Figure 9: RCS variation for the 3D finite-length circular PEC cylinder illuminated by an incident plane wave with $\varphi$${}^{inc}$ = 0${}^o$ and (a) $\theta$${}^{inc}$ = 90${}^o$ (xy-plane), (b) $\theta$${}^{inc}$ = 90${}^o$ (xz-plane), and (c) $\theta$${}^{inc}$ = 75${}^o$ (xz-plane), calculated via the reference staircase FDTD solution with $\mathrm{\Delta }$ = $\lambda$/72, the staircase NS-FDTD method only (without PI treatment, $\mathrm{\Delta }$ = $\lambda$/12, $\lambda$ = 1 m), the PI scheme ($\mathrm{\Delta }$ = $\lambda$/12) of [23], and the proposed PI model ($\mathrm{\Delta }$ = $\lambda$/12).

Figure 10: Antenna impedance, Z${}_{a\mathit{}n\mathit{}t}$ = R + jX, versus the wire length L for (a) $\mathrm{\Delta }$ = L/21 and a = 0.001 m, (b) $\mathrm{\Delta }$ = L/41 and a = 0.001 m, (c) $\mathrm{\Delta }$ = L/21 and a = 0.002 m, and (d) $\mathrm{\Delta }$ = L/41 and a = 0.002 m, computed through the proposed PI model and compared to the results from the NEC2 software and the staircase FDTD method.

2) Study of a thin-wire antenna impedance

The modified PI model for a thin-wire antenna of length L is elaborately described in Appendix B. The antenna is placed along the z-axis at the center of a 10$\mathrm{\Delta }$-wide PI region, surrounded by a NS-FDTD area terminated by a 20$\mathrm{\Delta }$-thick PML. To calculate the input impedance of the antenna, $Z_{\mathrm{ant}}=R+\mathrm{j}X$, we use (B11) and (B12) from Appendix B, while, in our simulations, we consider a normalized wavelength (i.e., $\lambda$ = 1 m); hence, the cell size is set to $\mathrm{\Delta }$ = L/21 or L/41. Moreover, for the combined PI and NS-FDTD area, $\mathrm{\Delta }$t is the same as its FDTD counterpart. The antenna is driven by a center feed with a $\mathrm{\Delta }$-gap on the z-axis and two wire radii are examined, namely a = 0.001 m and 0.002 m. Lastly, our reference solution is obtained via the NEC2 software [27], based on the method of moments (MoM) [28], as well as the staircase FDTD technique [16] with a properly fine lattice. As detected from the antenna impedance results of Fig. 10, the new PI model agrees better with the NEC2 than the FDTD method. This proves that the suggested formulation enables the use of the thin wire model [16], [25], [26] with the NS-FDTD algorithm, a remarkable advantage for the analysis of 3D curved configurations with demanding structural details.

V. CONCLUSION

The PI model for Maxwell’s equations is suitable to treat various arbitrarily-shaped objects when combined with FDTD schemes. To this objective, an efficient PI algorithm, equivalent to the NS-FDTD concept, has been presented in this paper for 2D and 3D problems. To keep this equivalence, the new PI model uses a pair of suitably tailored basic and complementary integral paths. So, its numerical accuracy and stability are the same as those of its NS-FDTD counterparts. Adapting and incorporating the featured PI model in the NS-FDTD analysis of realistic configurations, it has been shown that it offers high precision results for non-orthogonally-shaped structures, even on coarse grids. Thus, its combination with the NS-FDTD method increases the applicability of the latter and leads to notable computational savings.

APPENDIX A

The 3D NS-FDTD operator $d_x^{(0)}$ is expressed as

$$d_x^{(0)}={\alpha }_1{d}_x^{(1)}+{\alpha }_2{d}_x^{(2)}+{\alpha }_3{d}_x^{(3)},$$

(A1)

with

$$d_x^{(1)}f(x,y,z)=f(x+\mathrm{\Delta }/2,y,z)-f(x-\mathrm{\Delta }/2,y,z),$$

(A2)

$d_x^{(2)}$	$f(x,y,z)$	(A3)
$=$	${\displaystyle \frac{1}{4}}[f(x+\mathrm{\Delta }/2,y+\mathrm{\Delta },z+\mathrm{\Delta })-f(x-\mathrm{\Delta }/2,y+\mathrm{\Delta },z+\mathrm{\Delta })$
		$+f(x+\mathrm{\Delta }/2,y+\mathrm{\Delta }z-\mathrm{\Delta })-f(x-\mathrm{\Delta }/2,y+\mathrm{\Delta }z-\mathrm{\Delta })$
		$+f(x+\mathrm{\Delta }/2,y-\mathrm{\Delta },z+\mathrm{\Delta })-f(x-\mathrm{\Delta }/2,y-\mathrm{\Delta },z+\mathrm{\Delta })$
		$+f(x+\mathrm{\Delta }/2,y-\mathrm{\Delta },z-\mathrm{\Delta })-f(x-\mathrm{\Delta }/2,y-\mathrm{\Delta },z-\mathrm{\Delta })],$

		$d_x^{(3)}f(x,y,z)$	(A4)
		$={\displaystyle \frac{1}{4}}[f(x+\mathrm{\Delta }/2,y+\mathrm{\Delta },z)-f(x-\mathrm{\Delta }/2,y+\mathrm{\Delta },z)$
		$+f(x+\mathrm{\Delta }/2,y-\mathrm{\Delta },z)-f(x-\mathrm{\Delta }/2,y-\mathrm{\Delta },z)$
		$+f(x+\mathrm{\Delta }/2,y,z+\mathrm{\Delta })-f(x-\mathrm{\Delta }/2,y,z+\mathrm{\Delta })$
		$+f(x+\mathrm{\Delta }/2,y,z-\mathrm{\Delta })-f(x-\mathrm{\Delta }/2,y,z-\mathrm{\Delta })],$

where $\alpha$${}_1$ + $\alpha$${}_2$ + $\alpha$${}_3$ = 1. Note that operators $d_{y,z}^{(0)}$ can also be similarly derived. Parameters $\alpha$${}_{1,2,3}$ are given by

$${\alpha }_1={\eta }_1+{\eta }_2/3+{\eta }_3/2\text{and}{\alpha }_2={\eta }_2/3\text{,}$$

(A5)

for

$${\eta }_1=7/15-17{(k\mathrm{\Delta })}^2/864\text{and}{\eta }_2={\eta }_1-2{\gamma }_0+1\text{,}$$

(A6)

on condition that $\eta$${}_1$ + $\eta$${}_2$ + $\eta$${}_3$ = 1. A detailed description for the derivation of operators $d_{x,y,z}^{(0)}$, parameters $\alpha$${}_{1,2,3}$ and coefficient $\gamma$${}_0$ can be found in [1–3].

APPENDIX B

The modified integral paths for a thin-wire antenna are shown in Fig. B1. Thus, to compute the ${\partial }_t{\int }_S\mathbf{H}\cdot d\mathbf{S}$ quantity, the integral areas are $S=S_B={\mathrm{\Delta }}^2-\pi a^2/4$ and $S$ $=S_c=2{\mathrm{\Delta }}^2-\pi a^2/2$, as given in Fig. B1 (a). Moreover, after defining unit vector $\mathbf{t}$, we can write

Figure B1: Modified integral paths for a straight thin wire of radius a [black solid line: basic path; blue dashed line: complementary path]. In particular, integral path (a) near the wire at the xy-plane, (b) for a wire end on the z-axis, (c) near the wire along the z-axis (only the basic path is used), and (d) for a feed point with a $\mathrm{\Delta }$-gap on the z-axis. For simplicity, coordinate z in (a) and y in (b)-(d) are omitted, while the same procedure can be utilized at the yz-plane.

$$\mathbf{t}=\left(-{\mathbf{e}}_x+{\mathbf{e}}_y\right)/\sqrt{2},$$

(B1)

		${\tilde{E}}_y^{\prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2)$	(B2)
		$=\left[E_x(x,y-\mathrm{\Delta }/2){\mathbf{e}}_x+E_y(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }){\mathbf{e}}_y\right]\cdot \mathbf{t},$

		${\tilde{E}}_y^{\prime \prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2)$	(B3)
		$=\left[E_x(x-\mathrm{\Delta },y-\mathrm{\Delta }/2){\mathbf{e}}_x+E_y(x-\mathrm{\Delta }/2,y){\mathbf{e}}_y\right]\cdot \mathbf{t}.$

Since, E${}_x$(x, y – $\mathrm{\Delta }$/2) $\propto$ 1/x and E${}_y$(x – $\mathrm{\Delta }$/2, y) $\propto$ 1/y,

		${\displaystyle {\int }_{\mathit{\text{a⃝}}}^{\mathit{\text{d⃝}}}}\tilde{\mathbf{E}}\cdot d\mathbf{l}={\displaystyle \frac{\sqrt{2}\mathrm{\Delta }}{2}}{\displaystyle {\int }_a^{\sqrt{2}\mathrm{\Delta }/2}}{\tilde{E}}_y^{\prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2)l^{-1}dl$	(B4)
		$+{\displaystyle \frac{\sqrt{2}\mathrm{\Delta }}{2}}{\displaystyle {\int }_a^{\sqrt{2}\mathrm{\Delta }/2}}{\tilde{E}}_y^{\prime \prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2)l^{-1}dl$
		$={\displaystyle \frac{\mathrm{\Delta }}{\sqrt{2}}}\ln \left({\displaystyle \frac{\mathrm{\Delta }a}{\sqrt{2}}}\right)[{\tilde{E}}_y^{\prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2)$
		$+{\tilde{E}}_y^{\prime \prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2)].$

For the complementary path of Figs. B1 (b) and (d),

$${\tilde{E}}_z^{\prime }(x-\mathrm{\Delta }/2,z-\mathrm{\Delta }/2)=-E_x(x,z-\mathrm{\Delta }/2)\cos (\pi /4),$$

(B5)

$${\tilde{E}}_x^{\prime }(x-\mathrm{\Delta }/2,z+\mathrm{\Delta }/2)=E_x(x,z+\mathrm{\Delta }/2)\cos (\pi /4),$$

(B6)

where E${}_x$(x, y $\pm$ $\mathrm{\Delta }$/2) $\propto$ 1/x apart from E${}_x$(x, y + $\mathrm{\Delta }$/2) in Fig. B1 (b). In this manner, one obtains

$$\begin{array}{cc}\hfill {\int }_{\mathit{\text{a⃝}}}^{\text{4⃝}}{\tilde{\mathbf{E}}}^{\prime }\cdot d\mathbf{l}& =\frac{\sqrt{2}\mathrm{\Delta }}{4}{\int }_{\sqrt{2}a}^{\sqrt{2}\mathrm{\Delta }/2}{\tilde{E}}_y^{\prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2)l^{-1}dl\hfill \\ & =\frac{\sqrt{2}\mathrm{\Delta }}{4}\ln \left(\frac{\mathrm{\Delta }}{2a}\right){\tilde{E}}_y^{\prime }(x-\mathrm{\Delta }/2,y-\mathrm{\Delta }/2),\hfill \end{array}$$

(B7)

$$\begin{array}{cc}\hfill {\int }_{\text{3⃝}}^{\mathit{\text{c⃝}}}{\tilde{\mathbf{E}}}^{\prime }\cdot d\mathbf{l}& =\frac{\sqrt{2}\mathrm{\Delta }}{4}{\int }_{\sqrt{2}a}^{\sqrt{2}\mathrm{\Delta }/2}{\tilde{E}}_x^{\prime }(x-\mathrm{\Delta }/2,y+\mathrm{\Delta }/2)l^{-1}dl\hfill \\ & =\frac{\sqrt{2}\mathrm{\Delta }}{4}\ln \left(\frac{\mathrm{\Delta }}{2a}\right){\tilde{E}}_y^{\prime }(x-\mathrm{\Delta }/2,y+\mathrm{\Delta }/2).\hfill \end{array}$$

(B8)

On the other hand, for the basic path in Fig. B1(a),

${\displaystyle {\int }_{\text{4⃝}}^{\text{1⃝}}}\mathbf{E}\cdot d\mathbf{l}$	$={\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle {\int }_a^{\mathrm{\Delta }}}{E}_x(x,y-\mathrm{\Delta }/2)l^{-1}dl$	(B9)
		$={\displaystyle \frac{\mathrm{\Delta }}{2}}\ln \left({\displaystyle \frac{\mathrm{\Delta }}{a}}\right)E_x(x,y-\mathrm{\Delta }/2),$

which is similar to the 4⃝-3⃝ path in Fig. B1 (a), the 4⃝-1⃝ path in Figs. B1 (b)- B1 (d) and the 3⃝-2⃝ path in Figs. B1 (c)- B1 (d). So, the surface integral for the ${\partial }_t{\tilde{H}}_y$ term in Fig. B1 (c), using the basic path, is

$${\partial }_t{\int }_S{\overline{H}}_{y}dS\Rightarrow {\partial }_t{\overline{H}}_y(x,z)\frac{{\mathrm{\Delta }}^2}{2}\ln \left(\frac{\mathrm{\Delta }}{2a}\right),$$

(B10)

assuming that $H_y(x,z)\propto 1/x$. Note that all integral calculations, except for (B1)-(B10), treat field values on the path as a constant. Then, to evaluate the ${\partial }_t{H}_{y,z}$ derivatives, we can use (10), after replacing $\mathrm{\Delta }$ with $s_k(\mathrm{\Delta })$ through $\mu {\partial }_t{\int }_S\mathbf{H}\cdot d\mathbf{S}=-{\int }_C\mathbf{E}\cdot d\mathbf{l}$, while for the ${\partial }_t\mathbf{E}$ term, the basic path is employed. The feed point with the ${\mathrm{\Delta }}^{-}$ gap along $E_z(x-\mathrm{\Delta }/2,z)$ in Fig. B1 (d), is incorporated as

$${\int }_S\left(\epsilon {\partial }_t\mathbf{E}+\sigma \mathbf{E}+{\mathbf{j}}_{\text{source}}\right)\cdot d\mathbf{S}={\int }_C\mathbf{H}\cdot d\mathbf{l},$$

(B11)

where $\sigma$ is the electric conductivity and j${}_{s\mathbf{}o\mathbf{}u\mathbf{}r\mathbf{}c\mathbf{}e}$ the driving current. Thus, referring to Fig. $\mathrm{B}$1 (d) [16], [25], [26], the antenna impedance at the feed point is calculated by

$$Z_{\text{ant}}=R+\mathrm{j}X=\frac{-{\int }_{\text{4⃝}}^{\text{3⃝}}\mathbf{E}\cdot d\mathbf{l}}{{\int }_C\mathbf{H}\cdot d\mathbf{l}},$$

(B12)

with R the resistance and X the reactance.

REFERENCES

[1] J. B. Cole, “A high accuracy FDTD algorithm to solve microwave propagation and scattering problems on a coarse grid,” IEEE Trans. Microw. Theory Tech., vol. 43, no. 9, pp. 2053–2058, 1995.

[2] J. B. Cole, “A high-accuracy realization of the Yee algorithm using non-standard finite differences,” IEEE Trans. Microw. Theory Tech., vol. 45, no. 6, pp. 991–996, 1997.

[3] J. B. Cole, “High-accuracy Yee algorithm based on nonstandard finite differences: New developments and verifications,” IEEE Trans. Antennas Propag., vol. 50, no. 9, pp. 1185–1191, 2002.

[4] J. B. Cole, “High-accuracy FDTD solution of the absorbing wave equation, and conducting Maxwell’s equations based on a nonstandard finite-difference model,” IEEE Trans. Antennas Propag., vol. 52, no. 3, pp. 725–729, 2004.

[5] K. Taguchi, T. Ohtani, T. Kashiwa, and Y. Kanai, “Characteristics of evanescent waves in the nonstandard FDTD method,” IEEE Trans. Magn., vol. 43, no. 4, pp. 1313–1316, 2007.

[6] T. Ohtani, K. Taguchi, T. Kashiwa, Y. Kanai, and J. B. Cole, “Nonstandard FDTD method for wideband analysis,” IEEE Trans. Antennas Propag., vol. 57, no. 8, pp. 2386–2396, 2009.

[7] T. Ohtani and Y. Kanai, “Coefficients of finite difference operator for rectangular cell NS-FDTD method,” IEEE Trans. Antennas Propag., vol. 59, no. 1, pp. 206–213, 2011.

[8] N. Okada and J. B. Cole, “Nonstandard finite difference time domain algorithm for Berenger’s perfectly matched layer,” Applied Computational Electromagnetics Society (ACES) Journal, vol. 26, no. 2, pp. 153–159, 2011.

[9] T. Ohtani and Y. Kanai, “Characteristics of boundary model in the 2-D NS-FDTD method,” IEEE Trans. Magn., vol. 48, no. 2, pp. 191–194,2012.

[10] J. B. Cole and N. Okada, “High accuracy models for source terms in the nonstandard FDTD algorithm,” in Proc. Int. Symp. Electromagn. Theory, art. no. 24PM2C-01, pp. 1098–1100, 2013.

[11] T. Ohtani, Y. Kanai, and N. V. Kantartzis, “A 4-D subgrid scheme for the NS-FDTD technique using the CNS-FDTD algorithm with the Shepard method and a Gaussian smoothing filter,” IEEE Trans. Magn., vol. 51, no. 3, art. no. 7201004, 2015.

[12] T. Ohtani and Y. Kanai, “An enhanced Total-field/Scattered-field scheme for the 3-D nonstandard finite-difference time-domain method,” IEEE Trans. Magn., vol. 52, no. 3, art. no. 7204705,2016.

[13] J. B. Cole and S. Banerjee, Computing the Flow of Light: Nonstandard FDTD Methodologies for Photonics Designs. Bellingam, WA: SPIE Press,2017.

[14] J. Jose, S. K. Simon, J. Kizhakooden, J. Andrews, and V. P. Joseph, “Nonstandard FDTD realization of radiation behaviour of epsilon negative metamaterial corner reflector antenna,” in Proc. IEEE 13${}^{t\mathit{}h}$ Int. Congr. Artif. Mat. Novel Wave Phenom.-Metamat., pp. 178–180, Rome, Italy, 2019.

[15] J. B. Cole, R. Katouf, and S. Banerjee, “Nonstandard finite difference time domain methodology for harmonics generation in nonlinear dielectrics,” in Proc. 2021 Int. Appl. Comput. Electromagn. Soc. Symp., Online, 2021.

[16] A. Taflove and S. Hagness, Computational Electrodynamics: The Finite- Difference Time-Domain Method. 3rd Edition, Chs. 7, 8, 10, and 14. Norwood, MA: Artech House, 2005.

[17] A. Taflove, Advances in Computational Electrodynamics: The Finite-Difference Time-Domain Method. Norwood, MA: Artech House, 1998.

[18] D. M. Sullivan, Electromagnetic Simulation Using the FDTD Method. New York: IEEE Press, 2000.

[19] K. S. Kunz and R. J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics. Ch. 13, New York: CRC Press, 1993.

[20] C. J. Railton, I. J. Craddock, and J. B. Schneider, “The analysis of general two-dimensional PEC structures using a modified CPFDTD algorithm,” IEEE Trans. Microw. Theory Tech., vol. 44, no. 10, pp. 1728–1733, 1996.

[21] S. Dey and R. Mittra, “A locally conformal finite-difference time-domain (FDTD) algorithm for modeling three-dimensional perfectly conducting objects,” IEEE Microw. Guided Wave Lett., vol. 7, no. 9, pp. 273–275, 1997.

[22] T. Ohtani, Y. Kanai, and N. V. Kantartzis, “A rigorous path integral scheme for the two-dimensional nonstandard finite-difference time-domain method,” in Proc. COMPUMAG 2019, art. no. PA-M4–9, Paris, France, 2019.

[23] T. Ohtani, Y. Kanai, and N. V. Kantartzis, “A nonstandard path integral model for curved surface analysis,” Energies, vol. 15, art. no. 4322, 2022.

[24] T. Ohtani, Y. Kanai, and N. V. Kantartzis, “An integral representation model for the nonstandard finite-difference time-domain scheme,” in Proc. COMPUMAG 2021, art. no. OD1–4, Online,2022.

[25] S. Watanabe and M. Taki, “An improved FDTD model for the feeding gap of a thin-wire antenna,” IEEE Microwave Guided Wave Lett, vol. 8, no. 4, pp. 152–154, 1998.

[26] R. M. Makinen, J. S. Juntunen, and M. A. Kivikoski, “An improved thin-wire model for FDTD,” IEEE Trans. Microw. Theory Tech., vol. 50, no. 5, pp. 1245–1255, 2002.

[27] NEC: Numerical Electromagnetics Code, Available online: https://www.nec2.org/

[28] J. J. H. Wang, Generalized Moment Methods in Electromagnetics: Formulation and Computer Solution of Integral Equations. John Wiley & Sons: USA, 1991.

BIOGRAPHIES

Tadao Ohtani received the B.S. and M.S. degrees in electrical and electronic engineering from Toyohashi University of Technology, Japan, in 1983 and 1985, respectively, and received Ph.D. degree in electrical and electronic engineering from Kitami Institute of Technology, Japan, in 2005. From 1985 to 2011, he worked as a researcher at Nagoya Aerospace Systems of Mitsubishi Heavy Industries, Ltd. Currently, he is an independent researcher. His research interests include numerical analysis of the electromagnetic scattering fields for aircraft design via the FDTD and the NS-FDTD method.

Yasushi Kanai (Fellow, ACES) received the B.S., M.S. degree in engineering, and Ph.D. degree in information engineering from Niigata University, Japan, in 1982, 1984, and 1989, respectively.

He worked as a research engineer at Alps Electric Co., Ltd., from 1984 to 1992, where he developed magnetic recording heads via numerical methods. In 1992–1995, he was an associate professor at Department of Information Engineering, Niigata University. In 1995, he joined the Engineering Department, Niigata Institute of Technology, Kashiwazaki, Japan, where he is currently a professor. In 2002–2003, he was at Florida International University, Miami, FL, as a visiting scholar. He has authored/co-authored more than 190 peer-reviewed journal papers, more than 270 international conference records, more than 260 national conference records, and several book chapters. He specializes in micromagnetic analysis both in energy-assisted magnetic recording heads and media and in wave propagation via NS-FDTD analysis.

Nikolaos V. Kantartzis received the Diploma and Ph.D. degrees in electrical and computer engineering from the Aristotle University of Thessaloniki, Thessaloniki, Greece, in 1994 and 1999, respectively. In 1999, he joined the Department of Electrical and Computer Engineering, Aristotle University of Thessaloniki, where he is currently a professor. He has authored/coauthored 4 books, more than 190 peer-reviewed journal papers, and more than 300 publications in conference proceedings. His main research interests include computational electromagnetics, EMC, metamaterials, graphene, antenna design, and waveguide systems.

ACES JOURNAL, Vol. 39, No. 3, 189–200
doi: 10.13052/2024.ACES.J.390304
© 2024 River Publishers