—This work discusses the Finite-Difference Time-Domain (FDTD) technique to simulate an electromagnetic wave assuming one, two and three dimensions. The propagation medium is assumed to be a free space bounded by two absorbing boundaries, perfect matched layer (PML) and perfect electric conductor (PEC). The FDTD-1D is considered in free space while FDTD-2D and 3D are considered both in free space and in a free spae-medium consisting of dielectric sphere and cylinder in the center. In this case, we model the incident and the scattered electromagnetic ﬁelds reﬂected back from hitting the dielectric cylinder and sphere. Moreover, the simulation starts by generating an electromagnetic pulse either in the middle or at one end of the medium and this pulse can be either Gaussian or sinusoidal. For the FDTD- 3D, an antenna dipole is assumed to be the source generator of the electromagnetic pulse. We also provide the analytic solutions to conﬁrm the accuracy of the FDTD technique. Index Terms —FDTD, PEC, PML, Scattered ﬁelds. I. I NTRODUCTION A. A BRIEF HISTORY OF THE FDTD METHOD The Finite-Difference Time-Domain (FDTD) method was origi- nally proposed by Kane S.Yee in the seminal paper he published in 1966 [1]–[17]. Yee proposed a discrete solution to Maxwells equations based on central difference approximations of the spatial and temporal derivatives of the curl-equations [18]–[25]. The novelty of Yees approach was the staggering of the electric and magnetic ﬁelds in both space and time in order to obtain second-order accuracy. Yee derived a full three-dimensional formulation, and he validated the method with two-dimensional problems. Yees method went mostly un-noticed for nearly a decade.Finally,in 1975, Taﬂove and Brodwin applied Yees method to simulate the scattering by dielectric cylinders and biological heating, and in 1977, Holland applied it to predict the currents induced on an aircraft by an electromagnetic pulse (EMP). Holland summarized the reason for the delay of the application of Yees method: The algorithm was ﬁrst described by Yee about ten years ago. At that time, computers did not exist which could implement this scheme practically in three dimensions. This is no longer the case. It is now possible to mesh a problem space of interest into a 303030 grid and have random access to the resulting 162,000 ﬁeld quantities. These advances in computer technology have prompted Longmire to credit Yees algorithm as being the fastest numerical way of solving Maxwells equations. The growth of the FDTD method since the late 1970s can certainly be paralleled to the ongoing advances in computer technology. In fact, with todays massively parallel computers, a grid dimension of 3:0003:0003:000, which requires nearly 1.5 terabytes of memory, is feasible. Longmires quote above summarizes another reason for its growth. The FDTD method is a fast algorithm. If Nis the total number of degrees of freedom in the three-dimensional space, each time-iteration only requires O(N)ﬂoating-point operations. The caveat is that the discrete mesh has to ﬁll the full three-dimensional space. Thus, the number of degrees of freedom grows cubically with the linear dimension of the problem domain. Another reason for the wide-spread use of the FDTD method is Elyes Balti is with the Wireless Networking and Communications Group, Department of Electrical and Computer Engineering, The University of Texas at Austin, Austin, TX 78712 USA e-mail: ebalti@utexas.edu.its simplicity. One is reminded of Chef Auguste Gusteaus signature statement in Pixars movie Ratatouille: Anyone can cook!. Similarly, one can say: Anyone can FDTD!. The FDTD method is a wonderfully simple method that can be taught at the undergraduate or early graduate level. Yet, it is capable of solving extremely sophisticated engineering problems. The maturation of the FDTD method can be attributed to a number of pioneers who have made a host of contributions to advance Yee’s base FDTD algorithm. Each innovation has played a key role in advancing the FDTD method to where it is today. One set of key contributions revolve around the ability of the FDTD method to simulate unbounded problems. To model such problems, the discrete domain must be truncated via a reﬂectionless absorbing boundary.In 1981,G.Mur proposed a second-order accurate absorbing boundary condition (ABC) that helped to resolve this issue. Later, Higdon introduced a class of boundary operators that was more versatile. Betz-and Mittra extended Higdon’s ABC to absorb evanescent waves. While such ABC’s were widely used, the dynamic range of the FDTD method as well as the range of applications was still limited by this class of absorbing boundaries. A much more accurate absorbing boundary based on a Perfectly Matched Layer (PML) absorbing medium was proposed by J.-P. Berenger. In fact, the reﬂection error of the PML can be 3 or 4 orders of magnitude smaller than that offered by an ABC. The PML was also capable of truncating unbounded media that was lossy, inhomogeneous, dispersive, non- linear, or anisotropic. Further developments of Berenger’s PML also rendered it capable of absorbing evanescent waves and near ﬁelds.The tradeoff of using the PML is an extension of the mesh and additional degrees of freedom in the PML region. Thus, it is more computationally intensive than a local ABC. However, in the mid 1990’s, these resources were available even with commodity computers. The direct result of this innovation was that the FDTD method with PML absorbing boundaries was capable of being applied to a much broader range of problems. One of the advantages of the FDTD solution of Maxwells equations is that arbitrary media types can inherently be modeled. The FDTD method naturally accommodates inhomogeneous and lossy media. More complex media types, such as frequency-dependent dispersive, anisotropic, bi-anisotropic, chiral, or non-linear media can also be accommodated within the FDTD method. Signiﬁcant effort has been made by researchers to develop accurate and efﬁcient algorithms to model such media within the context of the FDTD method. An excellent summary of handling complex media within the FDTD method is provided by F. L. Teixeira. As the FDTD method has matured, the number of application areas it has impacted has diversiﬁed. Initially, the FDTD method was applied primarily to classical areas in electromagnetics, including electromagnetic scattering, electromagnetic compatibility, antennas, microwave circuits, and wave propagation. More recently, it has also been applied a diversity of other areas, including biomedical engineering, electromagnetic environmental hazards, ground penetrat- ing radar, photonics, biophotonics, plasmonics, photovoltaics, nano- optical storage devices, and seismic detection. The FDTD method is continuing to develop and expand. Current thrust areas include highorder FDTD methods,more general gridding techniques, uncon- ditionally stable schemes, and multiphysics applications. With these 2 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION endeavors and many others in progress, as well as continued advances in computing. B. LIMITATIONS OF THE FDTD METHOD While uncovering a historical perspective of the FDTD method, the previous section presented various strengths of the FDTD method. Of course, the method is not without its weaknesses, and it is certainly necessary to underline what these are. That is,it is important to understand the limitations of the FDTD method in order to determine when the method is most appropriately applied. One weakness of the FDTD method is that it requires a full discretization of the electric and magnetic ﬁelds throughout the entire volume domain. There are many instances when the FDTD method is forced to model a signiﬁcant amount of white-space, i.e., open space with no inhomogeneities. An example could be the electromagnetic scattering of a cluster of perfectly conducting spheres. The white space would be the volume inside each sphere, as well as the region in-between the spheres and the region separating the spheres and the absorbing boundary. The wider the separation between the spheres, the larger the percentage of white space. Another example is a printed circuit in a layered media where the traces being modeled constitute a very small percentage of the volume. The fact that the FDTD method is fully explicit can also be a weakness. While this can be considered a strength since a linear system of equations is not required to be solved, it can be a weakness if the time-step becomes inordinately small. The time- step is limited by the smallest geometric feature in the model. As a consequence, models with ﬁne, electrically small geometric features can have very small time-steps. This can result in a large number of time-iterations. This poses challenges when applying the FDTD method to large problems with ﬁne geometric features that must be modeled with high ﬁdelity. The FDTD method provides a broad-band simulation. As a conse- quence, a broad frequency response can be analyzed via a single simulation. Again, in many circumstances, this can be considered a strength. However, there are instances when only a narrow band response is desired. As a consequence, a frequency domain simula- tion may be more efﬁcient. Furthermore, when modeling complex materials, if the material constants are only known over a narrow frequency band, this characterization must be approximated over a broadband within the FDTD method. This could be as simple as the loss tangent of a substrate, or the effective material properties of a metamaterial. The FDTD method can also be challenged when the system under test has a very high Q. Consequently, the time-domain simulation may take a very long time to reach a steady state. Often, this is due to narrow band resonances that slowly decay (if at all). Thus, very long simulation times can result. For some problems, this can be mitigated by using methods such as the Generalized Pencil of Functions to either quantify the resonances, or to extrapolate the signal. In other cases, a frequency-domain simulation about the resonant frequencies may be more expedient. The FDTD method is also restricted by orthogonal gridding.This can be overcome by the use of subcell modeling techniques as well as sub-gridding and non-uniform grid methods. However, the local ﬁelds can still lose accuracy. Also, with methods such as subcell modeling, there is a level of uncertainty as to the exact shape of the local boundary. Thus, it is important to understand the level of accuracy of the near ﬁelds that can be achieved by subcell models for the geometry at hand. Higher-order FDTD methods have not shared the success of higher-order algorithms posed using other methodologies (such as discontinuous Galerkin methods). Most higher-order FDTD methods posed to date are based on an extended stencil namely,they are based on high-order difference approximations that employ points that span multiple grid cells. This makes it difﬁcult to model simple geometries that involve jump discontinuities in the materials. These methods are also challenged by structures with ﬁne geometric detail. II. Y EEALGORITHM The FDTD algorithm as ﬁrst proposed by Kane Yee in 1966 employs second-order central differences. The algorithm can be summarized as follows: 1) Replace all the derivatives in Amperes and Faradays laws with ﬁnite differences. Discretize space and time so that the electric and magnetic ﬁelds are staggered in both space and time. 2) Solve the resulting difference equations to obtain ’update equations’ that express the (unknown) future ﬁelds in terms of (known) past ﬁelds. 3) Evaluate the magnetic ﬁelds one time-step into the future so they are now known (effectively they become past ﬁelds). 4) Evaluate the electric ﬁelds one time-step into the future so they are now known (effectively they become past ﬁelds). 5) Repeat the previous two steps until the ﬁelds have been obtained over the desired duration. III. O NE-DIMENSIONAL SIMULATION TABLE I: Electromagnetic Parameters Symbol Q UANTITY B magnetic ﬂux density H magnetic ﬁeld strength D electric ﬁeld density E electric ﬁeld strength 0 permittivity of free space 0 permeability of free space r relative permittivity r relative permeability  conductivity  intrinsic impedance c speed of light A. FDTD Formulation Considering Maxwell’s curl equations in free space as follows: @E @t=1 0r H (1) @H @t=1 0r E (2) Given that along this work one dimension is assumed, the components ExandHyare only considered to facilitate the calculus and the other components are set to zeros. Given that the wave propagates in zdirection, the curl equations are reduced as follows: @Ex @t=1 0@Hy @z(3) @Hy @t=1 0@Ex @z(4) IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 3 To discretize these equations in time and space, we should apply the central difference method for temporal and spatial derivatives. The above aligns can be written as follows: En+1 2x(k)En1 2x(k) t=1 0Hn y k+1 2 Hn y k1 2 x(5) Hn+1 y k+1 2 Hn y k+1 2 t=1 0En+1 2x(k+ 1)En+1 2x(k) x (6) where tandxare the time and the space step, respectively. nand kare the number of the time and space steps, respectively. An instant tis given byt=nt, while a distance zis given byz=kx. The discretized equations can be formulated into an iterative algorithm as follows: En+1 2x(k) =En1 2x(k)t 0x Hn y k+1 2 Hn y k1 2 (7) Hn+1 y k+1 2 =Hn y k+1 2 t 0x  En+1 2x(k+ 1)En+1 2x(k) (8) Note thatExandHydiffer by a large order of magnitude because of0and0. To ﬁx this problem, we will roughly normalize the magnitude of the electric ﬁeld so that it becomes comparable to the amplitude of the magnetic ﬁeld. The new expression of the electric ﬁeld can be obtained by [ ?, eq. (1.5)]: eE=r0 0E (9) This new expression of the electric ﬁeld will be used for the analytic solution in order to get a result in agreement with the FDTD method. As we mentioned earlier, the space is equally dicretized into cells of the same size. For a given cell of size x, the time step is given by: t=x 2c(10) wherecis the speed of light. In this work, we create a Gaussian pulse in the middle and once the pulse is launched at t0, the wave propagates in the two directions (positive and negative zdirections) until it hits the boundaries, then it reﬂects back. Given that the boundaries are absorbing, the waves will not be reﬂected back. The ﬁrst boundary can be modelled as follows: En x(0) =En2 x(1) (11) The electric ﬁeld can be given at the second boundary as follows: En x(KE) =En2 x(KE1) (12) where KE is the length of the propagation medium. B. Analytic Solution Consider a Gaussian pulse as follows: f(t) = exp(t2) (13) The electromagnetic ﬁeld is assumed to propagate in both positive and negative zdirections. If the time pulse is Gaussian, the spectral density of this pulse is also Gaussian. The temporal Fourier transform of the initial pulse is given by: F(!) =r  exp !2 4 (14)Note that the pulse in frequency domain should be multiplied by a carrier or a propagator factor eikzto move in the positive zdirection and a carrier e+ikzto move in the negative zdirection. To get the expression of the electric ﬁeld in time and space (positive zdirection), the next step is to take the inverse of the Fourier transform of the spectral pulse multiplied by the carrier as follows: E+(z;t) =1 21Z 0F(!)eikze+i!td! (15) The electric ﬁeld E(z;t)propagating in the negative zdirection can be derived by replacing the propagating factor in eq. (15) by e+ikz. Finally, the total electric ﬁeld in time and space is given by: E(z;t) =E+(z;t) +E(z;t) (16) A Gaussian pulse is created at the instant t0in the middle of the space problem taking the following form: f(T) = exp" 0:5t0T spread2# (17) The Fourrier transform of the above pulse is given by: F(!) =spreade(!spread )2 2ei!t0(18) The electric ﬁeld in the positive z direction can be given by: E+(z;t) =spread2e(t0t)2 2spread 2eikz(19) The electric ﬁeld in the negative z direction can be given by: E(z;t) =spread2e(t0t)2 2spread 2e+ikz(20) Figures 1 and 2 show clearly the agreement between the FDTD and the analytic solution for different time steps and space. As we see, the wave starts from the middle and then it propagates in the two directions until the instant T= 250. At the instant T= 150, the wave moves to other positions and it is still far from the boundaries. At the instant T= 225, the wave starts hitting the two boundaries and it roughly overtakes the two ends at the instant T= 250 and then it totally vanishes from the space problem because it is absorbed by the two boundaries. IV. T WO-DIMENTIONAL SIMULATION A. Free Space Medium with Absorbing Boundaries We will start with the normalized Maxwells equations as follows: @eD @t=1p00reH (21) eD(!) = r(!)eE(!) (22) @H @t=1p00reE (23) For the TM mode, the above equations reduced as follows: @Dz @t=1p00@Hy @x@Hx @y (24) D(!) = r(!)E(!) (25) @Hx @t=1p00@Ez @y(26) 4 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 0 20 40 60 80 100 120 140 160 180 200 FDTD cells-1-0.500.51Electric Field Ex (V/m)T = 100 0 20 40 60 80 100 120 140 160 180 200 FDTD cells-1-0.500.51Magnetic Field Hy (A/m)T = 100 Fig. 1: FDTD-1D simulation of a Gaussian pulse in free space after 100 time steps. The pulse originated in the center and travels outward (KE = 200, t0= 40, spread = 12). 0 20 40 60 80 100 120 140 160 180 200 FDTD cells-0.500.511.52Electric Field Ex (V/m)T = 0 0 20 40 60 80 100 120 140 160 180 200 FDTD cells-0.500.511.5Electric Field Ex (V/m)T = 150 0 20 40 60 80 100 120 140 160 180 200 FDTD cells-0.500.511.5Electric Field Ex (V/m)T = 225 0 20 40 60 80 100 120 140 160 180 200 FDTD cells-0.500.511.5Electric Field Ex (V/m)T = 250 Fig. 2: Simulation of an FDTD program with absorbing boundary conditions. Notice that the pulse is absorbed at the edges without reﬂecting anything back. @Hy @t=1p00@Ez @x(27) The discretization of the above equations leads to: Dn+1=2 z (i;j)Dn1=2 z (i;j) t=1p00  Hn x i+1 2;j Hn y i1 2;j x! 1p00  Hn x i;j+1 2 Hn x i;j1 2 x!(28)Hn+1 x i+1 2;j Hn x i+1 2;j t=1p00 En+1=2 z (i;j+ 1)En+1=2 z (i;j) x(29) Hn+1 y i+1 2;j Hn y i+1 2;j t=1p00 En+1=2 z (i+ 1;j)En+1=2 z (i;j) x(30) Fig. 3 shows clearly the propagation of a Gaussian pulse generated in the middle of the space problem. We observe how the wave IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 5 0 600.20.4 60Electric Field Ez (V/m) 400.6 50T = 20 0.8 401 30 2020 100 0 0 600.20.4 60Electric Field Ez (V/m) 400.6 50T = 30 0.8 401 30 2020 100 0 0 600.20.4 60Electric Field Ez (V/m) 400.6 50T = 40 0.8 401 30 2020 100 0 0 600.20.4 60Electric Field Ez (V/m) 400.6 50T = 50 0.8 401 30 2020 100 0 Fig. 3: FDTD-2D simulation of a Gaussian pulse in free space medium with absorbing boundaries in various time instants. propagates in all directions till it hits the boundaries and then totally vanishes. B. The Perfectly Matched Layer (PML) According to Sullivan’s book, the PML parameters can be pre- sented as follows: Fig. 4: Parameters related to the perfectly matched layer (PML). Maxwell’s equations can be updated as follows: Dn+1=2 z (i;j) =gi3(i)gj3(j)Dn1=2 z (i;j) + 0:5gi2(i)gj2(j)  Hn y i+1 2;j Hn y i1 2;j Hn y i;j+1 2 +Hn y i;j1 2 (31)wheregi2;gi3;gj2;gj3are given by: gi2 =gj2 =1 1 +t 20 gi3 =gj3 =1t 20 1 +t 20 The curl of the electric ﬁeld can be given by: rE=En+1=2 z (i+ 1;j)En+1=2 z (i;j) (32) The incident magnetic ﬁeld in the Y-direction can be given by: In+1=2 Hy i+1 2;j =In1=2 Hy i+1 2;j +rE (33) The total magnetic ﬁeld in the Y-direction can be evaluated recur- sively as follows: Hn+1=2 y i+1 2;j =gi3 i+1 2 Hn y i+1 2;j 0:5gi2 i+1 2 rEgj1(j)In+1=2 Hy i+1 2;j (34) wheregi1is given by: gi1 =t 20 The incident magnetic ﬁeld in the X-direction is given by: In+1=2 Hx i;j+1 2 =In1=2 Hx i;j+1 2 +rE (35) 6 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION The total magnetic ﬁeld in the X-direction can be obtained by: Hn+1 x i;j+1 2 =gj3 j+1 2 Hn x i;j+1 2 + 0:5gj2 j+1 2 rE+gj1(j)In+1=2 Hy i;j+1 2 (36) Fig. 5 shows an FDTD-2D simulation of a sinusoidal source is initiated at a point that is offset ﬁve cells from the center of the problem space in each direction. As the wave reaches the perfectly matched layer (PML), which is eight cells on every side, it is absorbed. The effectiveness of the PML is apparent in the bottom ﬁgure because the contours would not be concentric circles if the outgoing wave was partially reﬂected. C. Total/Scattered Field Formulation In order to simulate a plane wave in a two-dimensional FDTD program, the problem space will be divided into two regions, the total ﬁeld and the scattered ﬁeld (in below ﬁgure). The two primary reasons for doing this are (i) the propagating plane wave should not interact with the ABCs and (ii) the load on the ABCs should be minimized. These boundary conditions are not perfect, that is, a certain portion of the impinging wave is reﬂected back into the problem space. By subtracting the incident ﬁeld, the amount of the radiating ﬁeld hitting the boundary is minimized, thereby reducing the amount of error. Fig. 6: Total ﬁeld/scattered ﬁeld of the two-dimensional space. Fig. 7: Every point is in either the total ﬁeld or the scattered ﬁeld. There are three places that must be modiﬁed: TheDzvalue at j = ja or j = jb: Dz(i;ja) =Dz(i;ja) + 0:5Hx;inc ja1 2 (37) Dz(i;jb) =Dz(i;jb)0:5Hx;inc ja1 2 (38)TheHxﬁeld just outside at j = ja or j = jb: Hx i;ja1 2 =Hx i;ja1 2 + 0:5Ez;inc(ja) (39) Hx i;jb+1 2 =Hx i;jb+1 2 0:5Ez;inc(jb) (40) Hyjust outside at i = ia and i = ib: Hy ia1 2;j =Hy ia1 2;j 0:5Ez;inc(j) (41) Hy ib1 2;j =Hy ib1 2;j + 0:5Ez;inc(j) (42) D. A Plane Wave Impinging on a Dielectric Cylinder To simulate a plane wave interacting with an object, we must specify the object according to its electromagnetic propertiesthe dielectric constant and the conductivity. For instance, suppose we are simulating a plane wave striking a dielectric cylinder 20 cm in diameter, which has a dielectric constant speciﬁed by the parameter epsilon and a conductivity speciﬁed by the parameter sigma. Fig. 9: Simulation of a plane wave striking a dielectric cylinder. The ﬁelds scattered from the cylinder are the only ﬁelds to leave the total ﬁeld and strike the PML. The simulation of a plane wave pulse hitting a dielectric cylinder with r= 30 and= 0.3 is shown in Fig. 10. After 25 time steps, the plane wave has started from the side; after 50 time steps, the pulse is interacting with the cylinder. Some of the pulse passes through the cylinder, and some of it goes around it. After 100 steps, the main part of the propagating pulse is being subtracted from the end of the total ﬁeld. E. Analytic Solution For the case of an electromagnetic wave propagating in free space medium bounded by PML, the propagator factor now consists of kx andky(2 directions). It can be given by: Space Propagator = exp [i(kxx+kyy)] (43) After reproducing the analysis of FDTD-1D for both XandY directions, the analytic solution can be derived. Regarding the case of a free space medium equipped with a dielectric cylinder, the scattered ﬁeld can be derived as follows: ETM s(r) = (E0 i^vi)k ki;1X m=1Jm(ki;;a) H(2) m(ki;a)(j)mexp(jm i) NH m;ki;z(; ;z ) (44) IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 7 -0.5 60 600Electric Field Ez (V/m) 4050T = 40 400.5 30 2020 100 0 -0.5 60 600Electric Field Ez (V/m) 4050T = 100 400.5 30 2020 100 0 Fig. 5: FDTD-2D simulation of a sinusoidal pulse with PML in various instants. -0.5 600 60Electric Field Ez (V/m) 40500.5T = 30 401 30 2020 100 0 -0.5 600 60Electric Field Ez (V/m) 40500.5T = 60 401 30 2020 100 0 -0.5 600 60Electric Field Ez (V/m) 40500.5T = 90 401 30 2020 100 0 -0.5 600 60Electric Field Ez (V/m) 40500.5T = 115 401 30 2020 100 0 Fig. 8: Simulation of a plane wave pulse propagating in free space. The incident pulse is generated at one end and subtracted at the other end. HTM s(r) =j 1X m=1mJm(ki;) H(2) m(ki;)MH m;ki;z(; ;z ) (45) wherem,MH m;ki;z,NH m;ki;z, m;ki;zare given by: m=(E0 i^vi)k ki;(j)mexp(jm i) MH m;ki;z=1 kr[ H m;ki;zz] NH m;ki;z=1 kr[MH m;ki;z] H m;ki;z(; ;z ) = [AH(1) m(kr) +BH(2) m(kr)] exp[jm jkzz]WhereH(1) mandH(2) mare the Hankel functions of the ﬁrst and second kind, respectively. V. T HREE -DIMENSIONAL SIMULATION A. Free Space Simulation Starting with Maxwell’s equations, we have: @eD @t=1p00reH (46) eD(!) = r(!)eE(!) (47) 8 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION -0.5 600 60Electric Field Ez (V/m) 40500.5T = 25 401 30 2020 100 0 -0.5 600 60Electric Field Ez (V/m) 40500.5T = 50 401 30 2020 100 0 -0.5 600 60Electric Field Ez (V/m) 40500.5T = 75 401 30 2020 100 0 -0.5 600 60Electric Field Ez (V/m) 40500.5T = 100 401 30 2020 100 0 Fig. 10: Simulation of a plane wave impinging on a dielectric cylinder. The cylinder is 20 cm in diameter and has a relative dielectric constant of 30 and a conductivity of 0.3 S/m. @H @t=1p00reE (48) After expanding the above equations, we have: @Dx @t=1p00@Hz @y@Hy @z (49) @Dy @t=1p00@Hx @z@Hz @x (50) @Dz @t=1p00@Hy @x@Hx @y (51) @Hx @t=1p00@Ey @z@Ez @y (52) @Hy @t=1p00@Ez @x@Ex @z (53) @Hz @t=1p00@Ex @y@Ey @x (54) After taking the ﬁnite-difference approximations on the equations (51) and (54), we have: Dn+1=2 z i;j;k +1 2 =Dn+1=2 z i;j;k +1 2 +t xp00  Hn y i+1 2;j;k+1 2 Hn y i1 2;j;k+1 2 Hn+1 x i+1 2;j+1 2;k +Hn+1 x i+1 2;j1 2;k (55)Hn z i+1 2;j+1 2;k =Hn z i+1 2;j+1 2;k +t xp00  En+1=2 x i+1 2;j+ 1;k En+1=2 x i+1 2;j;k En+1=2 y i+ 1;j+1 2;k +En+1=2 y i;j+1 2;k (56) Fig. 11 shows the propagation of the Ezfrom the dipole in the XY plane level with the gap of the dipole. Of course, there is radiation in the Z direction as well. This illustrates a major problem in three- dimensional simulations: unless one has unusually good graphics, visualizing three dimensions can be difﬁcult. B. A Plane Wave Impinging on a Dielectric Sphere in Free Space Medium Bounded by PML The Maxwell’s equations in this case can be updated as follows: rH=Hn y i+1 2;j;k+1 2 Hn y i1 2;j;k+1 2 Hn+1 x i+1 2;j+1 2;k +Hn+1 x i+1 2;j1 2;k(57) The incident electric ﬁeld density is given by: In Dz i;j;k +1 2 =In1 Dz i;j;k +1 2 +rH (58) Dn+1=2 z i;j;k +1 2 =gi3(i)gj3(j)Dn+1=2 z i;j;k +1 2 +0:5gi2(i)gj2(j) rH+gk1(k)In+1=2 Dz i;j;k +1 2(59) IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 9 0 400.05 30 400.1Electric Field Ez (V/m)T = 30 300.15 200.2 201010 0 0 0 400.05 30 400.1Electric Field Ez (V/m)T = 35 300.15 200.2 201010 0 0 0 400.05 30 400.1Electric Field Ez (V/m)T = 45 300.15 200.2 201010 0 0 0 400.05 30 400.1Electric Field Ez (V/m)T = 60 300.15 200.2 201010 0 0 Fig. 11:Ezﬁeld radiation from a dipole antenna in a three-dimensional FDTD program. Now that we have a program that generates a plane wave in three dimensions, we will want to start putting objects in the problem space to see how the plane wave interacts with them. In two dimensions, we chose a cylinder because we had an analytic solution with which we could check the accuracy of our calculation via a Bessel function expansion. It turns out that the interaction of a plane wave with a dielectric sphere can be determined by an expansion of the modiﬁed Bessel functions. Fig. 12: Total/scattered ﬁelds in three dimensions. The total scattered electric ﬁeld density in Y direction can be obtained at the two ends by: Dy i;j+1 2;ka =Dy i;j+1 2;ka 0:5Hx;inc(j)(60)Dy i;j+1 2;kb+ 1 =Dy i;j+1 2;kb+ 1 0:5Hx;inc(j) (61) The analytic solutions are the same as the case of FDTD-2D. VI. C ONCLUSION In this paper, we present a global framework analysis of FDTD 1D, 2D and 3D considering various combination of the medium properties such as absorbing boundaries, perfect matched layer, dielectric cylinder, etc. We provide pictorial presentations of the wave propagation in various instant to conﬁrm the behavior of the wave in the medium with the analytic solution. In addition, we present the design of the patch antenna and the optimal parameters that provide the required antenna properties. Capitalizing on these parameters, we derive the antenna ﬁgure of merits such as the pattern, the directivity at both azimuthal and elevation planes, the input impedance, the VSWR, the return loss and the jS11j. As an extension of this work, we intend to design an array antenna with powerful bandwidth around 60 GHz to introduce it the world of the mmWave communication. In addition, to increase the spectral efﬁciency, we should consider a massive MIMO that consists of latge number of small antenna (large densiﬁed array) and the design of such structure is very challenging which will be the topic of our next research work.