# Electric field Discontinuity-Considered Effective-Permittivities and Integration-Tensors for the three-dimensional Finite-Difference Time-Domain method.

1. INTRODUCTION

A popular numerical tool for solving electromagnetic problems is the Finite-Difference Time-Domain (FDTD) method. FDTD discretizes the problem space using cells with particular values of electric-permittivity and magnetic-permeability and calculates Maxwell's equations in a leap-frog manner. Inherent problem with this approach is that non-planar interfaces between two different materials can only be approximated. One of the crudest forms of approximations is the staircase method because of the jagged nature of the approximated interface between two materials.

There have been various efforts to overcome the errors from approximations. Kaneda et al. [2] proposed a method that obtains effective dielectric constants for cells containing two different materials. This method first calculates a sectional dielectric constant by the area-weighted sum of two different materials. Subsequently, the average of the inverse of this value is calculated along the direction normal to the section. Finally, the inverse of this average obtained in each axis is used as the effective dielectric constant. Yu and Mittra [3] proposed much simpler approach by only considering the length-weighted sum of the dielectric constants of two materials along each axis. Dey and Mittra [4] used the volume-averaged sum of the dielectric constants of two materials.

Above three methods do not take account for the interface normal directions. Recently, a new method called subpixel smoothing [5-7] was introduced that accounts for the normal direction of the interface between two materials as well as their volumes ratio within a cell. The contour-path method [1,8] which is the two-dimensional analogy of this paper also considers the interface normal direction and the volumes ratio within a cell. However, the derivations of the two methods are quite different. Subpixel smoothing starts by postulating a dielectric tensor that was selected assuming that parallel and normal electric field experience different dielectric constants [9]. This dielectric tensor is then smoothed where the smoothing kernel becomes a Dirac delta function when the kernel envelope normal to the interface reaches zero. This smoothed dielectric tensor is finally tested by substituting in the first-order corrected angular frequency which is obtained by expanding an eigensolution of the vector Helmholtz equation in terms of the small change in the dielectric so that the limiting value is well defined [10]. On the other hand, contour-path method uses the integral form of the Ampere and Faraday laws and obtains a tensor form of the dielectric constant in a straight forward manner. The obtained dielectric tensor is unique based on series of steps that involve the consideration of boundary conditions between dielectric discontinuities. Subpixel smoothing was successfully applied to three-dimensions but the contour-path method was only applied in two-dimensions. This paper presents a three-dimensional extension of the contour-path method and it is shown that the two methods agree in the effective dielectric tensor. This finding is important because it gives justification to the selection of the effective dielectric tensor in the subpixel smoothing method.

While paying much attention to the discontinuity of the electric field vector, we found it to be natural to apply similar procedures to the line integration of the electric field in the integral form of the Faraday's law. Previous contour-path method only considered the electric field discontinuity in the integral form of the Ampere's law. The gain in this second part was found to be significant in the numerical computation of the scattering problems investigated when the discretization error is large.

In this paper, we have compared our method to the staircase, volume-averaged and subpixel smoothing methods. We will be only considering dielectric, nonmagnetic materials.

2. CONTOUR-PATH METHOD FOR DIELECTRICS

This section gently introduces the derivation of contour-path method in three-dimensions. We will first discuss Ampere and Faraday laws in three-dimensions and then move on to effective-permittivities. We then present the electric field Discontinuity-Considered Effective-Permittivities and Integration-Tensors (DC-EP&IT) that are the main topics of this paper.

Figure 1 illustrates Yee's cell. A cell stores two vectors that represent an electric field and a magnetic field. Rather than storing the three values of a vector at the cell's center, Yee's scheme was to store the values at shifted locations. This shifted storage is very well suited for representing Ampere and Faraday laws. For nonmagnetic dielectric materials, Ampere and Faraday laws become

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

where E and H are electric and magnetic fields; [epsilon] and [mu] are electric-permittivity and magnetic-permeability; dS and dL are differential surface and length elements whose directions are normal to the surface and parallel to the curve, respectively.

[FIGURE 1 OMITTED]

For example, the Ampere law in the z-axis at z = (k + 1/2) [DELTA]h can be written as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

where [DELTA]h and [??] are grid distance and the unit vector in the z-axis, z-axis, respectively. Notice we have deferred the dot product after the integration. And also note that since the Yee's cell only stores one field component at each grid location, the missing components are obtained by simply interpolating four adjacent Yee's grid locations. This scheme is used throughout this paper.

When we calculate the terms for the Faraday law, we need to calculate the line integration of the electric field. For example in the y-axis, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

where [??] and [??] are the unit vectors in the x and z-axes, respectively. We have similarly deferred the dot product after the integration.

Notice we have preserved the line integration in Eq. (3) but did not do so in Eq. (2). The purpose is for improving the accuracy. Recall that we are dealing with nonmagnetic materials where magnetic-permeability is uniform and there is no discontinuity in the magnetic fields. Therefore, we can safely evaluate the line integration in Eq. (2). However, there is a dielectric discontinuity in Eq. (3) and doing so would increase errors.

In order to evaluate Eqs. (2) and (3), we need to understand an important distinction between the electric field value stored at a grid location and the actual electric field value at that particular spatial location. The important understanding is that they are different. Electric field value stored represents that of the whole cell. Mohammadi et al. [1,8] assumed they are equal. Initial attempt assuming they are equal resulted in large errors and a new interpretation based on volume-averaged approach was derived that is explained below.

2.2. Representative Electric Field Vector

The volume integral of the electric field in a given time step inside a cubic cell with length [DELTA]h centered at i, j, k + 1 / 2 is E[|.sup.n+1/2.sub.i,j,k+1/2] [([DELTA]h).sup.3]. When the cell is composed of two dielectric materials of constants [[epsilon].sub.1] and [[epsilon].sub.2], the electric field will be discontinuous as [E.sub.1] and [E.sub.2]. According to the boundary conditions, we know that [[epsilon].sub.1] x [E.sub.1] = [[epsilon].sub.2]nn x [E.sub.2] since [D.sub.[perpendicular to],1] = [D.sub.[perpendicular to],2] across dielectrics, where nn is a dyadic, formed by juxtaposing a pair of interface normal vectors n and n. Also, we can write (I - nn) x [E.sub.1] = (I - nn) x [E.sub.2] since [E.sub.[parallel],1] = [E.sub.[parallel],2] across dielectrics, where I is an identity dyadic. Thus with the two boundary conditions and [V.sub.1] and [V.sub.2] representing the fractional volumes of material 1 and 2, the volume integral can be represented either entirely with [E.sub.1] or [E.sub.2]. That is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

We have placed [perpendicular to],[parallel] below some terms to highlight perpendicular and parallel terms. The above can be simplified as follows depending on which form of the electric field is used.

(I + ([[epsilon].sub.1] / [[epsilon].sub.2] - 1)[V.sub.2]nn) x [E.sub.1][|.sup.n+1/2.sub.i,j,k+1/2] [([DELTA]h).sup.3], (5)

(I + ([[epsilon].sub.1] / [[epsilon].sub.2] - 1)[V.sub.1]nn) x [E.sub.2][|.sup.n+1/2.sub.i,j,k+1/2] [([DELTA]h).sup.3]. (6)

There is no significance in which form to use: Eq. (5) or Eq. (6), when dealing with planar incident waves and planar dielectric interfaces. However, when non-planar waves and interfaces are concerned, it becomes more appropriate to use Eq. (5) when [V.sub.1] > [V.sub.2] and Eq. (6), otherwise. The reason is that the normal and parallel terms are derived approximations obtained by multiplying with nn and thus it makes more sense to keep these less exact terms multiplied by the smaller volume.

In summary, Eqs. (5) and (6) equal to E[|.sup.n+1/2.sub.i,j,k+1/2] [([DELTA]h).sup.3], and they give the relations between what is actually stored at the grid versus [E.sub.1] and [E.sub.2].

2.3. Electric Field Discontinuity-considered Effective-permittivities

In FDTD, the electric field is stored at grid locations and the integration shown on the left of Ampere's law in Eq. (2) is calculated by the simple product of the electric field and the effective dielectric constant. In the conventional scheme, the effective dielectric constant is a single value. However, in order to consider the discontinuity in the electric field, this has to be a tensor. The effective dielectric tensor of the top (pink) cell illustrated in Figure 2 is [epsilon][|.sub.i,j,k+1/2] and those in the x and y axes are [epsilon][|.sub.i+1/2,j,k] and [epsilon][|.sub.i,j+1/2,k], respectively. Note that each cell simultaneously resides in two Yee's cells that are abutting in the z, x or y-directions. In summary, we will be using three dielectric tensors in a given Yee's cell. The staircase and volume-averaged methods implemented for comparisons also used these shifted permittivities in three axes, although they are single values not tensors.

There are several steps involved in deriving the Discontinuity-Considered Effective-Permittivities (DC-EP). We first apply the integration form of Ampere's law (See Appendix A1 for details) and integrate the result along the normal direction of the cross section. We repeat this process for x, y and z-axes. Subsequently, we combine the three scalar equations into one vector equation to get the following dielectric tensors, depending on which electric field (1 or 2) was used (See Appendix A2 for details). With some lengthy work, it can be shown that the two equations are actually identical and it does not matter which form to use.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

It is also interesting to see that the dielectric tensors in Eq. (7), although quite different as they stand, are identical to that of subpixel smoothing [5]. Similar derivations can be applied for the x and y components where the tensor forms are the same.

2.4. Electric Field Discontinuity-considered Integration-tensors

We can also apply similar technique to the line integration for the electric field crossing dielectric discontinuities. For example, the line integration of the electric field in the z direction can assume either

[DELTA]h((I + ([[epsilon].sub.1] / [[epsilon].sub.2] - 1) [L.sub.2]nn) x [((I + ([[epsilon].sub.1] / [[epsilon].sub.2] - 1) [V.sub.2]nn).sup.-1] x E[|.sup.n+1/2.sub.i,j,k+1/2]) x [??] (8)

or

[DELTA]h((I + ([[epsilon].sub.2] / [[epsilon].sub.1] - 1) [L.sub.1]nn) x [((I + ([[epsilon].sub.2] / [[epsilon].sub.1] - 1) [V.sub.1]nn).sup.-1] x E[|.sup.n+1/2.sub.i,j,k+1/2]) x [??],

depending which fraction of volume is greater. Derivations are given in Appendix A3. When [V.sub.1] > [V.sub.2], Eq. (8) is used and Eq. (9) is used otherwise. [L.sub.1] and [L.sub.2] represent the fractional lengths of material 1 and 2. Electric field Discontinuity-Considered Integration-Tensors (DC-IT) shown in Eqs. (8) and (9) are used to evaluate the right hand side integration terms in Faraday's law in Eq. (3). Similar derivations can be done for x and y components where the tensor forms are the same.

[FIGURE 2 OMITTED]

Since generally, [L.sub.1] [not equal to] [V.sub.1] and [L.sub.2] = [V.sub.2], Eq. (8) and Eq. (9) cannot be simplified to [DELTA]h [E.sub.z][|.sup.n+1/2.sub.i,j,k+1/2] and the equality we saw for the dielectric tensors in Eq. (7) does not exist between Eq. (8) and Eq. (9).

3. RESULTS

We have implemented the proposed DC-EP and DC-IT in a three-dimensional FDTD code. To amplify the calculation errors, we have chosen to simulate a scattering by 100 nm diameter spherical object with high (12/1) dielectric constant ratio (object to the medium) using very coarse to fine cell sizes experiencing a planar wave of wavelengths ranging from 400 nm to 1100 nm. Scattered-field formulation was used to introduce the incident beam to the simulated electromagnetic field. More information about our FDTD code can be found in Sung et al. [12].

Justification to the selected electric-permittivities is needed. It is well known that the wavelength of an electromagnetic wave compresses when dielectric constant increases. When uniform grid lengths are employed in the computational space, compression effectively makes the cell coarser as compared to that with lower dielectric constant. It would be reasonable to conclude that the effect of higher dielectric constant ratio would resemble that of coarser cells and it would also be gradual. Therefore, we believe a dramatic change in the errors based on a critical dielectric contrast ratio would be unlikely. The particular (12/1) permittivity ratio was chosen based on a high contrast example in a real application (a photonic crystal material [5]).

Accuracies of normal directions, fractional volumes and lengths are important in computing DC-EP and DC-IT. We have approximated them by sampling methods and details are described in Appendix B. It is true that the integration by using sampled points can possibly reduce the accuracy of the computation. We could make the integration analytic for specific surface equations but that would make the method much more difficult to implement because various singular cases and special cases need to be addressed. Note this is also a problem for subpixel smoothing where the dielectric and inversed dielectric constants as well as interface normals need to be computed. For the current results, we have used 19 sampling points on each axis. For general cases that would make only about 5% of maximum error on each axis and 15% of maximum error in 3-dimensions.

Based on high dielectric constant ratio and big cell size, discretization error was purposely made more prominent than other sources of errors such as phase velocity, scattering field formulation and the perfectly matched layers. For the results, we have used the direct electric field difference error,

Error = [absolute value of [E.sub.FDTD] - [E.sub.Mie]] / [absolute value of [E.sub.Mie]], (10)

averaged for all Yee's cell on a bounding box whose planes are at two grid distance away from the sphere along each axis. Analytical Mie calculations were done by a code derived from the MIEV package [13] where the spatial offset of each component of the electric field vector was considered. The log value of this error is plotted in Figure 3 for frequency at 1100 nm, as a function of cell size which equals to wavelength divided by [N.sub.Lambda]. It can be seen that the DC-EP&IT and subpixel smoothing show far better results than the staircase and volume-averaged methods. The staircase method is equivalent to the classical FDTD method. Not shown but similar tendencies were shown for [lambda] = 400 nm, [lambda] = 650 nm and [lambda] = 900 nm.

The FDTD method employs finite differences as approximations to both the spatial and temporal derivatives that appear in Maxwell's equations. An examination of a difference formula for the second derivative on an equidistant grid leads to O([([DELTA]h).sup.2]) errors, In other words, it would show perfect-quadratic convergence in ideal conditions. This implies that if [DELTA]h is reduced by a factor of 10, the error in the approximation should be reduced by a factor of 100. However, due to other sources of errors, quadratic convergence shown in Figure 3 is seldom achieved. In practical sense, most methods show perfect-linear O([DELTA]h) or slower convergence. DC-EP&IT and subpixel methods show quadratic convergence up to [N.sub.Lamda] = 40 and a linear convergence thereafter.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

The error difference between subpixel smoothing and DC-EP&IT is highlighted in Figure 4 at [lambda] = 400 nm, [lambda] = 650 nm, [lambda] = 900 nm and [lambda] = 1100 nm, where a positive value denotes DC-EP&IT is better in terms of accuracy. The error enhancement as compared to the subpixel method is quite large on average about 10% at low subdivisions and reaches as high as 15%. Not much gain at [lambda] = 400 nm is shown due to the relative fine cell size.

4. CONCLUSION

The discontinuity of the electric field in the dielectric discontinuity was carefully considered in the integration form of the Ampere and the Faraday laws to come up with the electric field Discontinuity-Considered Effective-Permittivities and Integration Tensors (DC-EP&IT). The accuracy of the proposed method in comparison with the staircase, volume-averaged and subpixel methods was evaluated for the Mie scattering problem. It showed that the DC-EP&IT performed better than the other methods especially when the discretization error was large. Especially, the benefit as compared to the subpixel method was clearly shown when the cell size was large. This makes DC-EP&IT very well suited for tackling large problems where the size of the cell cannot be made small due to the limit in the computer memory size.

APPENDIX A. ELECTRIC DISCONTINUITY-CONSIDERED DERIVATIONS

A.1. Sectional Effective-permittivities

To derive the sectional effective-permittivities in tensor forms we first start with an imaginary sectional cut of the Yee's cell at z = [z.sub.cut] as shown in Figure A1.

According to the boundary condition between two dielectric materials of constants [[epsilon].sub.1] and [[epsilon].sub.2], we know that [[epsilon].sub.1]nn x [E.sub.1] = [[epsilon].sub.2]nn x [E.sub.2] since [D.sub.[perpendicular to],1] = [D.sub.[perpendicular to],2] across dielectrics and (I - nn) x [E.sub.1] = (I - nn) x [E.sub.2] since [E.sub.[parallel],1] = [E.sub.[parallel]2] across dielectrics. Thus with [A.sub.1] and [A.sub.2] representing the fractional area of material 1 and 2, the left hand side surface integration in Eq. (2) in a difference form becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A2)

With [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A3)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A4)

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

[FIGURE A1 OMITTED]

In summary, depending which electric field value is used; Eq. (A2) and Eq. (A4) are the effective-permittivities, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], at the location z = [z.sub.cut].

Similar derivations can be obtained in x and y directions and the results are summarized in Eq. (A5). Notice each dielectric constant can have two forms depending which electric field is used. [??], [??] and [??] are unit vectors in x,y and z-axes, respectively.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A5)

[FIGURE A2 OMITTED]

A.2. Derivations of Discontinuity-considered Effective-permittivities

Equations (A2) and (A4) hold at arbitrary cut location for infinitesimally thin sheet. For example in the z-axis, we can have many of these within the cell centered at i, j, k + 1 / 2. Figure A2 illustrates one possible thin sheet located at the center of this cell. When Eq. (A2) form of effective dielectric tensor is used, these can be represented using an auxiliary variable w as

([DELTA]h / [DELTA]t [[epsilon].sub.1,z][|.sub.i,j,k+1/2+w] x ([E.sub.1][|.sup.n+1/2.sub.i,j,k+1/2+w] - [E.sub.1][|.sup.n-1/2.sub.i,j,k+1/2])) x [??] - [H.sub.x][|.sup.n.sub.i,j+1/2,k+1/2+w] + [H.sub.y][|.sup.n.sub.i,j- 1/2,k+1/2+w] + [H.sub.y][|.sup.n.sub.i+1/2,j,k+1/2+w] - [H.sub.y][|.sup.n.sub.i-1/2,j,k+1/2+w], (A6)

Where w [member of] [-1/2,+1/2].

If we assume ([E.sub.1][|.sup.n+1/2.sub.i,j,k+1/2+w] - [E.sub.1][|.sup.n-1/2.sub.i,j,k+1/2+w]) [approximately equal to] [E.sub.1][|.sup.n+1/2.sub.i,j,k+1/2] -[E.sub.1][|.sup.n-1/2.sub.i,j,k+1/2] is constant, by integrating along the direction z we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A7)

We can also do a similar integration along the x and y direction within the cell centered at i,j,k + 1/2 to get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A8)

"representative value" refers to an average value in the one dimensional integration. It becomes appropriate to discard the representative value notation if the cell is sufficiently small and we are considering the value at the mid-plane of the cell. It can be shown that the three integrations of the tensor dielectric constant in x, y and z directions are identical. Thus we can drop the x, y and z subscript and represent these integrations as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A9)

This is also true for the duel dielectric tensor

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A10)

Similarly by substituting Eqs. (6) and (A10) into (A12), we get

Using Eq, (A9), Eqs. and (A8) can be combined leading to a vector equation,

[DELTA]h / [DELTA]t [[epsilon].sub.1][|.sub.i,j,k+1/2] x ([E.sub.1][|.sup.n+1/2.sub.i,j,k+1/2] - [E.sub.1][|.sup.n- 1/2.sub.i,j,k+1/2]) = [H.sub.curl], (A11)

Where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

And similarly, we can show that for the duel

[DELTA]h / [DELTA]t [[epsilon].sub.2][|.sub.i,j,k+1/2] x ([E.sub.2][|.sup.n+1/2.sub.i,j,k+1/2] - [E.sub.2][|.sup.n- 1/2.sub.i,j,k+1/2]) = [H.sub.curl]. (A12)

Substituting Eqs. (5) and (A9) into (A11), we get

[DELTA]h / [DELTA]t (([V.sub.1][[epsilon].sub.1] + [V.sub.2][[epsilon].sub.2])I + [V.sub.2]([[epsilon].sub.1] - [[epsilon].sub.2])nn) x [(I + ([[epsilon].sub.1] / [[epsilon].sub.2] - 1)[V.sub.2]nn).sup.-1] x (E[|.sup.n+1/2.sub.i,j,k+1/2] - E[|.sup.n-1/2.sub.i,j,k+1/2]) = [H.sub.curl]. (A13)

Similarly by substituting Eqs. (6)and (A10) into (A12), we get

[DELTA]h / [DELTA]t (([V.sub.1][[epsilon].sub.1] + [V.sub.2][[epsilon].sub.2])I + [V.sub.1]([[epsilon].sub.2] - [[epsilon].sub.1])nn) x [(I + ([[epsilon].sub.2] / [[epsilon].sub.1] - 1)[V.sub.1]nn).sup.-1] x (E[|.sup.n+1/2.sub.i,j,k+1/2] - E[|.sup.n-1/2.sub.i,j,k+1/2]) = [H.sub.curl]. (A14)

With some lengthy work, we can show that the resulting combined dielectric tensors in Eqs. (A13) and (A14) are equal,

A.3. Derivations of Discontinuity-considered Integration-tensors

We can also apply similar technique to the line integration for the electric field crossing dielectric discontinuities as we did for the surface integration. For example, the line integration of the electric field in the z direction can assume either

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A15)

or

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A16)

depending which fraction of volume is greater. When [V.sub.1] > [V.sub.2], Eq. (A15) is used and Eq. (A16) is used otherwise. [L.sub.1] and [L.sub.2] represent the fractional lengths of material 1 and 2. Figure A3 illustrates two cases depending which material segment is longer. It is important to understand that the selection of the proper equation is not dictated by which segment is longer. Similar derivations can be achieved for other line integrations in x and y directions where the tensor forms are the same.

[FIGURE A3 OMITTED]

APPENDIX B. COMPUTING NORMALS AND AREAS

Exact computation of the normal directions used in Eqs. (7), (8) and (9) can be quite demanding. We show how we have done this by approximations. We first calculate the normals at the boundaries by intersecting the edges of integrating voxel with the dielectric surface and subsequently average these to obtain the normal. Intersection between a sphere and a line segment was done by using the code from Cychosz et al. [11]. Figure B1(a) shows a case when the dielectric surface crosses two oppositely located edges at one face of a voxel while Figure B1(b) shows that of two adjacent edges. In Figure B1(a), [bar.n] = ([n.sub.30] + [n.sub.12]) / 2 and in Figure B1(b), [bar.n] = ([n.sub.30] + [n.sub.23]) / 2. Besides these edges, other crossing edges within the voxel need to be considered, which is omitted here for brevity.

The fractional volumes in Eqs. (7), (8) and (9) are also numerically obtained by distributing uniformly sampled points in the integrating square and performing point location query with respect to the dielectric surface. Figure B2 illustrates the sampled points on regular grid points. The ratio of the number of inside points to the total number of sampled points multiplied by the area of the integrating square becomes [A.sub.2]. [A.sub.1] is simply obtained by subtracting [A.sub.2] from the area of the integrating square. Afterwards, the areas along the normal direction at various altitudes are similarly computed and added. The sample points are {(u[DELTA]h/(2m+1), v[DELTA]h/(2m+1), w[DELTA]h/(2m+1))|u,v,w [member of] 0, [+ or -]1, [+ or -]2,...,[+ or -]m}, when the surface normal is in the z direction. Sample points in other directions can be similarly obtained. For the results in this paper, we have used m = 9.

[FIGURE B1 OMITTED]

[FIGURE B2 OMITTED]

ACKNOWLEDGMENT

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (grant number 2010-0028190) and by a grant from the institute of Medical System Engineering (iMSE) in the GIST, Korea.

Received 3 June 2011, Accepted 27 June 2011, Scheduled 6 July 2011

REFERENCES

[1.] Mohammadi, A., H. Nadgaran, and M. Agio, "Contour-path effective permittivities for the two-dimensional finite-difference time-domain method," Optics Express, Vol. 13, 10367-10381, 2005.

[2.] Kaneda, N., B. Houshmand, and T. Itoh, "FDTD analysis of dielectric resonators with curved surfaces," IEEE Transactions on Microwave Theory and Techniques, Vol. 45, 1645-1649, 1997.

[3.] Yu, W. and R. Mittra, "A conformal finite difference time domain technique for modeling curved dielectric surfaces," IEEE Microwave and Wireless Components Letters, Vol. 11, 25-27, 2001.

[4.] Dey, S. and R. Mittra, "A conformal finite-difference time-domain technique for modeling cylindrical dielectric resonators," IEEE Transactions on Microwave Theory and Techniques, Vol. 47, 1737-1739, 1999.

[5.] Farjadpour, A., D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, "Improving accuracy by subpixel smoothing in the finite-difference time domain," Optics Letters, Vol. 31, 2972-2974, 2006.

[6.] Oskooi, F., C. Kottke, and S. G. Johnson, "Accurate finite-difference time-domain simulation of anisotropic media by subpixel smoothing," Optics Letters, Vol. 34, 2778-2780, 2009.

[7.] Deinega, A. and I. Valuev, "Subpixel smoothing for conductive and dispersive media in the finite-difference time-domain method," Optics Letters, Vol. 32, 3429-3431, 2007.

[8.] Mohammadi, A., T. Jalali, and M. Agio, "Dispersive contour-path algorithm for the two-dimensional finite-difference time-domai method," Optics Express, Vol. 16, 7397-7406, 2008.

[9.] Lee, J.-Y. and N.-H. Myung, "Locally tensor conformal FDTD method for modeling arbitrary dielectric surfaces," Microwave and Optical Technology Letters, Vol. 23, 245-249, 1999.

[10.] Meade, R. D., A. M. Rappe, K. D. Brommer, J. D. Joannopoulos, and O. L. Alerhand, "Accurate theoretical analysis of photonic band-gap materials," Physical Review B, Vol. 48, 8434, 1993.

[11.] Cychosz, J. M. and W. N. Waggenspack, Jr., "Intersecting a ray with quadric surface," Graphics Gems III, 275-283, Academic Press Professional, Inc., 1992.

[12.] Sung, S.-Y. and Y.-G. Lee, "Trapping of a micro-bubble by non-paraxial Gaussian beam: Computation using the FDTD method," Optics Express, Vol. 16, 3463-3473, 2008.

[13.] Wiscombe, W. J., "Improved Mie scattering algorithms," Appl. Opt., Vol. 19, 1505-1509, 1980.

Y.-G. Lee *

Department of Mechatronics, Gwangju Institute of Science and Technology (GIST), 261 Cheomdan-gwagiro, Buk-gu, Gwangju 500712, Korea

* Corresponding author: Yong-Gu Lee (lygu@gist.ac.kr).
No portion of this article can be reproduced without the express written permission from the copyright holder.

Author: Printer friendly Cite/link Email Feedback Lee, Y.-G. Progress In Electromagnetics Research Report 9SOUT Aug 1, 2011 4849 Design and optimization of equal split broadband microstrip Wilkinson power divider using enhanced particle swarm optimization algorithm. Field synthesis in inhomogeneous media: joint control of polarization, uniformity and SAR in MRI [B.sub.1]-field. Dynamical systems Electric fields Tensors (Mathematics) Time-domain analysis