A Graphical User Interface for Scattering Analysis of Electromagnetic Waves Incident on Planar Layered Media

This paper introduces a MATLAB-based Graphical User Interface (GUI) which could help electromagnetics engineers and researchers who are interested in designing layered media for various applications. The paper begins with presenting the analysis method the program employs, continues by encountering specific considerations and techniques of implementation, and ends with providing different numerical examples. These examples show good efficiency of the program for analysis of diverse problems.


Introduction
Layered or composite Media have drawn lots of attention due to their wide spectrum of applications. They could serve as microwave and optical filters [1,2], coatings for cloaking and radar cross section (RCS) reduction [3][4][5][6], boosted transparency surfaces [7], high-reflection coatings [8], multilayer circuits [9], etc. This fact has not only resulted in myriads of researches pertaining their usages, but also has inspired many novel analysis methods, some of which are: transfer matrix method (TMM) [10,11], classical numerical techniques such as finite difference time domain (FDTD) [12] and integral equations [13], transmission line equivalence [14], and recursive methods [15,16]. This paper uses a relatively simple and efficient method which is based on transfer matrix method. The mathematical formulation is elaborated in [6] and will be restated briefly here with stressing some practical considerations. The reason we choose the transfer matrix method is firstly for its generality to solve nearly any related homogeneous theoretical and practical structure and secondly, for its clear and well-developed algorithm to be converted into MATLAB code. The provided algorithm is applicable for lossless and lossy media, anisotropic as well as isotropic ones, dispersive media, and complex structures composed of metamaterials and conventional media. The only restriction is that each layer should be homogeneous in every direction. Of course in special cases which the inhomogeneity is aligned with the structure's normal vector (here z axis), the structure may be viewed as infinite number (in practice, large number) of homogeneous layers and analyzed using current program. The proposed algorithm employs matrices and so is suitable to be realized by MATLAB which has a matrixbased coding. The program receives constitutive parameters of structure layers and incident wave properties namely frequency, incidence direction and polarization and gives the scattered wave parameters. The incident and scattered waves are all assumed to be planar in this research. After describing mathematical formulation, it will be verified using a few significant numerical examples.

Formulation
The geometry of the structure and incident and scattered wave vectors are shown in figure 1. The wave travels from left halfspace to the right halfspace that are assumed to be air (or some other isotropic media). Also, the right halfspace may be perfect electric conductor (PEC). This case is especially important when the structure is supposed to be used as a radar absorbing material (RAM) or anti-reflection coating. Each homogeneous layer in its most general form can be identified with its permittivity and permeability tensors (which are frequency dependent in common) and its thickness. We do not impose any limitation for permittivity and permeability tensors' elements in order to keep the generality of the problem. However, mostly in the literature, these tensors are assumed diagonal. The physical limitation that this assumption impose is that the layer's principal axis should be aligned with the normal vector of the structure. Note that the propagation direction lies in the x-z plane in figure 1 and this will be applied throughout our formulation. However, this does not reduce the generality of the problem as far as our formulation contains the general form of permittivity and permeability tensors. That's because if the propagation be not in x-z plane, we can take it in x-z plane and instead, rotate the whole structure around the z axis in order to keep the relative angle untouched. To apply such rotation, tensors should be updated by conversions in which ε , μ , R z , and  are the relative permittivity, relative permeability, rotation matrix about z axis, and propagation angle, respectively.
To know the overall scattering behavior of the structure, we should learn how planar waves propagate in any medium possessing aforementioned parameters. Assuming all plane wave fields as harmonic functions of space and time as [17] which in combination give the wave equation as: 1 2 0 where 0 k is the wave number in free space. Equation (1) is a vector equation and could also be represented by matrix form as: where κ is a matrix that transforms any vector like E to k E  and its value is: while the value of x k is the same in the whole space due to phase matching at the boundaries and equals to sin( ) k is the wave number in the left halfspace and  denotes the incident wave propagation angle .
Equation (2) is the matrix wave equation and in order to yield nonzero solution for electric field, 1 2 0 k   κμ κ ε should be singular. This condition leads us to equation: (4) which is a fourth order equation in terms of z k and has four solutions in general. Therefore, there are four wave vectors each of which expresses a different propagation mode. The overall field inside the medium is the linear superposition of all modes' fields as: where m E and ˆm E are the electric field excitation coefficient and unit vector of the mode "m" respectively and is the wave vector of that mode. The excitation coefficients are determined using the boundary conditions. To calculate the unit vectors, we return to equation (2). This matrix equation is in fact a set of three linear equations for E elements. Each of these equations is equation of a plane in ( , , ) x y z E E E space and the three planes intersect in one line that passes from origin and on which the E solution lies. Thus, the electric unit vector is the same as the unit direction vector of the line and equals to the normalized cross product of planes' normal vectors. In sum, assuming to be the r'th row of the matrix, the unit electric field would be: (6) for any r and  r if   r r . The value of above equation does not depend on the choice of r and  r , since three planes intersect in a single line. Equation (6) should be worked out for each zm k . The rest of calculations is the ordinary transfer matrix method.
To apply transfer matrix method, we need to define characteristic matrix of each layer that contains all the information about that layer and relates the field values of layer's two sides. The n'th layer characteristic matrix is defined as The matrix given in equation (7) should be calculated for the left and right halfspaces as well as layers and 0 Β and 1 N  Β will stand for them. To do this calculation, we use 1 2 3 4 x z   , and 4ˆĉ os( ) sin( ) for the left halfspace and similar relations for the right halfspace.
Next, we define the scattering matrix which relates the incident and scattered waves' fields in such a way that: (9) and could be calculated as: To take advantage of these relations, we define reflection and transmission matrices as:  (11) and: 13 14 11 12 where ij S is the ij'th component of scattering matrix. If the structure is PEC-backed, there is no transmission ( 0  T ) and from equation (11) the reflection matrix turns out to be: Now, the scattered fields could be obtained through: Equations (14) and (15) are the ultimate solutions and are able to be put into use by a computer code to solve the layered media scattering problems.

Algorithm of implementation
The proposed formulation is algorithmic and ready to use by a computer. The algorithm requests calculations which are basic matrix algebra and are embedded as basic operations in MATLAB. The only time-consuming equation to solve by numerical solution that MATLAB employs is equation (4).
For simpler and faster solution of equation (4), we use another procedure as explained subsequently.
To solve equation (4) (16) and the solutions of   6   k d is plotted in Figure 3.
The solution is in agreement with original article. However, the transfer matrix method fails for large values of 0 k d due to round-off errors as enormous exponential terms arise. For these values, other methods like recursive methods should be used. Figure 3: The GUI layout and its response for the first example.
This example demonstrates the accuracy of the proposed algorithm for smaller k0d values and its unstable response for higher frequencies.

Dispersive DPS-DNG bilayer structure
One great advantage of our algorithm is that its solution is valid for any sign of permittivity and permeability elements without extra considerations. In fact, the solutions of (4) satisfy the criteria of choosing correct sign [4] when the medium has negative values in its ε or μ . Therefore, we are able to analyze any class of DPS (double positive), DNG (double negative), ENG (epsilon negative), and MNG (mu negative) materials with the same equations. Combining these media in one structure, yields several interesting cases like absolute transparency. Figure 4 shows the reflectance (reflected wave power normalized to incident power) against frequency from a pair of DPS and DNG isotropic layers with thicknesses 1 and 2 millimeters that acts as a backward notch filter. The first layer is assumed to be a non-magnetic material with relative permittivity of 2.2 and the second layer is a lossless metamaterial with dispersion relations of table 1. Dispersion parameters are given in the figure legend. Agreement with [5] proves the potential of the algorithm to analyze highly dispersive metamaterial structures.

Perfectly matched layer
Perfectly matched layer (PML) is an extremely lossy medium that is used to absorb radiation in boundaries of numerical electromagnetic problems (mostly when finite difference time domain is adopted) to make the simulation space finite. The most common type of these media is Gedney's PML [20] which is defined as a uniaxial anisotropic material with permittivity elements of 0 for zz elements to model the Gedney's PML. Gedney's PML is matched to the air for any incidence angle, so is able to make better absorption relative to ordinary isotropic absorber with 0 1    j     parameters-which is perfectly matched only for normal incidence-when the incidence angle is not zero. However, for the normal incidence, the two cases are identical for equal conductivities. Figure 5 compares the isotropic and anisotropic (Gedney's) absorbers' power for killing TE-polarized electromagnetic wave with frequencies between 0-10 GHz and at angles of 0-90 degrees when they are backed by PEC. The layer's thickness and conductivity are assumed to be 1 (cm) and 1 (S/m). While the two act excellent for normal incidence, the anisotropic case is more efficient for oblique waves (as it can be seen from its concavity).

Conclusions
First, mathematical formulation based on the transfer matrix method for planar electromagnetic wave scattering analysis of planar layered structures consisting homogeneous materials that can be lossy, anisotropic, dispersive, and metamaterial was presented. The mathematical formulation was optimized for being coded by MATLAB scripting language and for easy and standalone use of the code, it was packaged as a graphical user interface. After considering the algorithm, the program was validated by three examples, and the good agreement with reports in resources, proved the integrity of the program.