{'fraction_non_alphanumeric': 0.03322917638280011, 'fraction_numerical': 0.011029754687062774, 'mean_word_length': 4.926616915422885, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reliable Nonuniform Discretization Algorithm for Fast and Accurate Hybrid Mode Analysis of Multilayered Planar Transmission Lines

A flexible and reliable full-wave modal integral method is proposed to efficiently characterize planar transmission structures printed on multilayered isotropic/anisotropic substrates. Based on the mathematical concept of operators used in electromagnetism, it consists in determining the modal inner products obtained through the Galerkin’s procedure via a proper choice of trial functions with metallic edge effects. To this aim, a fast and accurate nonuniform discretization algorithm is introduced for the first time, while using a new process to accelerate the convergence with regard to the number of areas of such inner products, thus significantly reducing the required CPU-time for planar transmission lines analysis. To demonstrate the efficiency of the proposed numerical integral approach, a successful comparison was achieved through a close agreement with published data.


INTRODUCTION
Widely used in RF/microwave systems, microwave planar transmission lines exhibit a hybrid propagation behavior due to their inhomogeneous propagation medium, which make their characterization highly challenging, resulting in complex boundary conditions, particularly at the airdielectric interface. Consequently, with the occurrence of increasingly complex structures in modern RF/microwave devices, the associated Helmholtz's equation becomes more and more difficult to solve.
To this aim, several powerful modeling numerical methods have been proposed to provide solutions to such propagation issues. However, to the best of the authors' knowledge, there is no existing approach that globally satisfies the key challenges encountered in the analysis of multilayered isotropic/anisotropic planar circuits such as high accuracy, minimum CPU time, and reduced memory requirement. Besides, we have to add the versatility of the selected method so that it could be applied to a large panel of circuits and topologies. Developing an approach that efficiently combines the above features is therefore critical. Several modeling methods have been proposed to address this issue. However, most of them suffer from large memory requirement and CPU time, particularly in characterizing asymmetrical planar lines in multilayer configurations [1]. Therefore, many numerical techniques have been retained to tackle propagation issues in multilayered structures such as the finite difference time domain (FDTD) method [2][3][4], spectral domain approach (SDA) [5][6][7], or multimode equivalent network (MEN) technique [8,9], to name a few. However, they suffer from some severe limitations. Thus, the FDTD method requires significant numerical efforts due to the complexity of its algorithm, based on the discretization of differential equations. As for the SDA formulation, it is not straightforward for inhomogeneous or truncated substrates because of the difficulty to set up the spectral dyadic Green's matrix. Regarding the MEN formulation, to the best of the authors' knowledge, this technique was only applied to isotropic substrates. Therefore, to support the ever-increasing use of multilayered isotropic/anisotropic planar transmission lines in hybrid mode, we opted for the full-wave integral technique combined with the formalism of mathematical operators used in electromagnetism (EM) and the method of moments involving the Galerkin's procedure in the modal domain. The rationale behind such formalism is to describe, in a more convenient way, the boundary conditions of the electromagnetic fields at the interfaces [10][11][12]. In fact, this technique has been successfully used to analyze lossy transmission lines with finite metallization thickness [13], symmetrical/asymmetrical transmission lines including isotropic/anisotropic substrates [1,14], truncated substrates [15,16], and multilayered substrates [17,18]. Note that the proposed approach has been developed from the equivalent electrical circuit concept that allows converting the continuity relations of EM field expressions into current-voltage relationships [10,11], thus leading to the determination of a key parameter, i.e., the admittance/impedance operator [1,14]. Thanks to the relative simplicity of the considered equivalent circuit, there is no restriction on the number of layers, which, indeed, significantly extends this integral method to the analysis of generic symmetrical/asymmetrical planar transmission structures with an arbitrary number of layers.
To better evaluate the field distribution or the current densities at the discontinuities due to the presence of metallized obstacles or abrupt changes in the dielectric medium, it is more appropriate to consider the fictitious propagation in the transverse direction instead of the longitudinal one (TEn and TMn modes). The above technique allows obtaining the admittance/impedance operator related to the considered structure to which the Galerkin's procedure was applied, leading to a dispersion matrix, whose elements are constituted of inner products (of basis functions with trial functions) deduced from the integral calculation [17]. Setting the determinant of the above matrix to zero ensures the existence of non-trivial physical solutions from which the propagation characteristics of the different modes can be deduced.
However, due to constraints related to the heavy numerical computation involved, this modal method requires truncating the series as well as the size of the dispersion matrix, which, in turn, implies to decide on an optimal choice of trial functions in order to reach the desired solution within a tolerable error range. Therefore, the reliability of such method is based on the proper choice of trial functions with respect to the region to be characterized, either metallic or insulating, according to the type of considered planar transmission structures printed on multilayered isotropic/anisotropic substrates. In fact, this region is chosen to minimize the numerical processing while ascertaining a suitable compromise between accuracy, speediness and memory size.
Furthermore, when trial functions taking into account metallic edge effects are involved, determining the inner products in the modal domain is quite laborious numerically. So, most designers have circumvented this problem by selecting trial functions with no metallic edge effects due to their simplicity and ease of integral computation [13,16,[19][20][21] such as rooftop or sinusoidal trial functions. In fact, these latter functions allow avoiding spurious solutions [22] but at the risk of reduced accuracy.
To efficiently resolve these challenging issues in the modal domain, we considered two types of trial functions. The first, the sinusoidal trial functions, lightens the analytical calculations but with a slower convergence with respect to the dispersion matrix size. The second, the trial functions with metallic edge effects, favor a faster convergence but at the cost of more complicated mathematical processing. Therefore, to address this issue while accurately determining the inner products, a new approach is proposed. It is based on a numerical integration using a nonuniform mesh in which the discretization step should be very small near the conducting strip edges. The proposed nonuniform discretization algorithm involves trapezoidal or quadratic elementary areas, while introducing for the first time a technique that significantly accelerates the convergence with regard to the number of areas. It is worth to mention that this numerical approach is generic, applicable to the analysis of any 2D or 3D multilayered planar structure; it suffices to change only the trial functions and the integration boundaries according to the study region (i.e., isolating or conducting region) of the considered structure.
Note that applying this modal technique to symmetrical lines with finite metallization thickness will involve a dispersion matrix containing 4-times the number of sub-matrices compared to the one with zero-thickness, hence a matrix two-times larger [23]. Therefore, to avoid complex calculations (and the inherent numerical errors it will implicate), we assumed conducting strip(s) with zero-thickness. This hypothesis can be considered as reasonable knowing that the metallization thickness is assumed much smaller than the transverse dimensions. In fact, the focus has been made mainly on speeding-up the convergence of the numerical integral computation.
Furthermore, according to the uniformity of the transmission structures along the propagation axis, a 2D analysis was considered. We also assumed homogeneous, isotropic/anisotropic and lossless dielectric substrates as well as perfectly conductive metallization with negligible thickness.
To demonstrate the efficiency of the proposed numerical integral approach with adequate choice of trial functions, successful comparisons were achieved between the obtained numerical results and published data.

MODAL METHOD DEVELOPMENT VIA MATHEMATICAL OPERATORS
The development of such modeling method via the formalism concept of mathematical operators helps reinforcing the systematic character of this technique thus involving a clear way to solve the boundary conditions [11,13]. The objective of this formalism is the evaluation of the overall admittance/impedance operator to which the Galerkin's procedure will be applied. In fact, such key operator represents a link between the tangential fields and the current densities. It involves an equivalent circuit of the transverse structure that includes a non-zero virtual source represented by a set of trial functions defined on the considered metallized discontinuity plane of the proposed transmission planar structure. For this purpose, let us consider a generic shielded multilayer unilateral symmetrical/asymmetrical transmission planar structure with arbitrary uniaxial bi-anisotropic stratified media above and below the metallized interface on which several metallic strips have been arbitrarily deposited, as shown in Fig. 1(a). This figure shows the cross sectional view of the mentioned multilayered general structure with its equivalent circuit defined either in the insulating or metallic domain. This structure is constituted of M layers and two regions with regard to the metallized interface, the above with M 1 layers and the below with M 2 layers, with M = M 1 + M 2 . Note that the uniaxial permittivity and permeability tensors for the ith dielectric layer can be expressed in the form of diagonal matrices such as which can be noted asε ri = [ε rxi , ε ryi , ε rzi ] andμ ri = [µ rxi , µ ryi , µ rzi ], with ε rxi = ε rzi = ε ri⊥ , µ rxi = µ rzi = µ ri⊥ , ε ryi = ε ri// , and µ ryi = µ ri// . Here, ε ri⊥ (µ ri⊥ ) and ε ri// (µ ri// ) are the respective transverse and longitudinal relative permittivities (permeabilities) with regard to the optical axis along the oy direction. The modal technique was developed by considering a fictitious propagation in the oy transverse direction instead of the real oz longitudinal direction. As shown in Fig. 1(b) and Fig. 1(c), the transverse section of the structure was modeled as a set of cascaded transmission lines terminated by short-circuited loads, noted "sc". This leads to an equivalent circuit that includes a transverse electric field " ⃗ E" (or current density " ⃗ J") described by trial functions defined in the slot(s) or on the metallic strip(s) of the structure.Ŷ I andŶ II are the admittance operators seen by the metallized interface for the respective upper and lower region. The Kirchoff's laws allow writing ⃗ As mentioned in [1,14,17], theŶ andẐ operators can be described in the form of diagonal matrices obtained from their spectral representation, thusẐ =Ŷ −1 . Their evaluation requires the involvement of a complete set of orthogonal basis functions ({|f n ⟩} n=0...N ) and the calculation of the mode admittances Y I n and Y II n of the respective upper and lower regions at the metallized interface (with N the number of modes).   Note that such basis functions [17,23] should satisfy the boundary conditions imposed by the shielding. So, the admittance operator can be obtained as: The impedance operator can be then deduced from the admittance operator as: where "e" and "h" indicate the TMn and TEn modes along the transverse direction (oy), respectively. Here, Y n and Z n represent the respective mode admittance and mode impedance viewed at this interface. |f n ⟩ ⟨f n | represents the projection operator on the basis vectors {|f n ⟩} n=0...N , knowing that the product of vector "bra" (⟨f n |) with vector "ket" (|f m ⟩) represents their inner product deducted from the integral calculation involving only the "x" variable such as: The matrix representation of the admittance/impedance operator requires the determination of the mode admittance/impedance obtained from the recursive formulation [17,18].
The application of the Galerkin's technique, for the asymmetrical structure with an arbitrary number of slots and metallic strips, requires the decomposition of the unknown electric field defined in the slots or current density on the metallic strips in terms of trial functions pondered by their coefficients such as: In the slots: On the strips: where K indicates the number of trial functions per component along the ox and oz directions, and I max is the number of slots or metallic strips. (Φ pxi , Φ qzi ) or (Φ ′ pxi , Φ ′ qzi ) are the transverse trial function components along ox and oz directions represented in the ith slot or metallic strip, respectively.
For structures involving one slot or one metallic strip, like finlines or microstrip lines, respectively ( Fig. 2(a) and Fig. 2(b)), the transverse electric fields or current densities are represented as follows: Note that (Φ px , Φ qz ) are the transverse trial functions defined in the slot of the symmetrical/asymmetrical finline and (Φ ′ px , Φ ′ qz ) the transverse trial functions defined on the metallic strip of the symmetrical/asymmetrical microstrip line.
The transverse current densities or electric fields are zero in the slot or on the metallic strip, respectively. They can be determined using the admittance or impedance operator as: or The application of the Galerkin's technique to Eq. (10) or (11), while involving the respective admittance or impedance operator, requires the use of testing functions identical to trial functions expressed in the slot or on the metallic strip, which allows obtaining a homogeneous system of algebraic equations according to the insulating or conducting domain. In fact, for each specific domain, we can obtain the vector containing the weighting coefficients of trial functions as well as the dispersion matrix (of dimension 2K × 2K) consisting of 4 sub-matrices composed of mode admittance/impedance and inner products of trial functions with basis functions [17,18].
For asymmetrical transmission configurations involving two slots at the dielectric-air interface, such as coplanar structures ( Fig. 3(a)), or two metallic strips like coupled microstrip structures ( Fig. 3(b)), the transverse electric fields or current densities can be represented as follows:  Here (Φ px1 , Φ qz1 ) and (Φ px2 , Φ qz2 ) are the respective trial functions defined in the first and second slot of the asymmetrical coplanar lines.
) are the respective trial functions defined on the first and second metallic strip of the asymmetrical coupled microstrip lines, with respective widths d 1 and d 2 . Knowing that the transverse current densities or electric fields are zero in the two slots or on the two metallic strips, respectively, these quantities can be determined using the admittance or impedance operator and Eq. (12) as: or The Galerkin's technique involving trial functions expressed in the two slots or on the two metallic strips leads to a homogeneous system of equations [17,18]. It should be noted that, for the symmetrical case (i.e., d 1 = d 2 = d), only half of the coplanar (or coupled microstrip) structure was considered in the analysis. Therefore, the unknown transverse electric fields (or current densities) can be expressed on the base of trial functions Φ px and Φ qz (or Φ ′ px and Φ ′ qz ) defined only in the right slot (or on the right metallic strip). In addition, transverse electric fields (or current densities) can be determined in terms of trial functions by using the expressions in Eq. (12), thus, the Galerkin's technique requires the same expressions as in Eq. (13) or (14). Of course, trial functions of such coplanar or coupled microstrip lines do not have the same expressions as those of finline or microstrip structures.
As for the size of the resulted dispersion matrix, it will be of 4K × 4K, constituted of 16 submatrices for asymmetrical coupled microstrip [14] and coplanar lines [1], and of 2K × 2K, composed of 4 sub-matrices for the symmetrical case and for symmetrical/asymmetrical finline and microstrip structures [17]. Indeed, the size of the global dispersion matrix depends on the number of slots (or metallic strips), i.e., 2I max K × 2I max K for asymmetrical lines with 4I 2 max sub-matrices, where I max indicates the number of slots (or metallic strips) [17]. Thus, the non-trivial solution, obtained by setting the determinant of the dispersion matrix to zero, allows evaluating, at a given frequency, the propagation parameters of the dominant and higher order propagation modes.

CHOICE OF TRIAL FUNCTIONS
The accuracy of the numerical results, obtained through the Galerkin's technique by solving the homogeneous system, depends mainly on the proper choice on trial functions, which also allows speedingup the convergence with regard to the number of trial functions per component (K) and modes (N ) while ensuring physical solutions [24]. Indeed, such correct selected functions involve a good approximation of the expanded transverse electric fields (or current densities) in the slot(s) (or on the metallic strip(s)) of the considered planar structure. However, this choice must respect some criteria as mentioned in [17], such as boundary, proportionality, and metallic edge effect conditions. In addition, it should respect the parity condition for the symmetrical structure. In fact, due to the presence of the electric wall (or the magnetic wall) on the symmetry plane, the x-and z-components of such transverse trial functions must be respectively both even and odd (or odd and even) with respect to the middle of the symmetry plane of the considered line [25,26]. Furthermore, this parity condition imposes that the nth mode should be even, i.e., n = {0, 2, 4, . . .} (or odd, i.e., n = {1, 3, 5, . . .}) according to the electrical (or magnetic) wall on the symmetry plane. Nevertheless, because it is difficult to simultaneously satisfy all the above conditions, we have to optimize the choice of these functions to fulfill these criteria as much as possible. For this purpose, a set of sinusoidal type trial functions with and without metallic edge effects have been selected, knowing that trial functions taken into account edge effects are known by their tendency to become infinite at the vicinity of the metallic edge. It involves a reduction of the dispersion matrix size obtained from the previous homogeneous system. Tables 1 and 2 summarize some trial functions that take into account metallic edge effects for shielded symmetrical/asymmetrical lines.
Note also that the transverse trial functions of symmetrical finlines and microstrip lines can be deducted by setting S 1 = (a − w)/2 (see Fig. 2) and by choosing the origin in the middle of those structures with regard to the x-axis. As for symmetrical coplanar and coupled microstrip lines, trial   Fig. 3), where the origin of the x-axis is on the top of these planar lines. Note that these latter functions must also verify the parity conditions described by the even and odd modes.
It should be highlighted that the above expressions do not seem to satisfy the proportionality condition, involving thus the existence of spurious propagation parameters [22] (described as infinite solutions) instead of physical ones. However, these non-physical solutions could be easily eliminated by choosing only the solutions obtained by the modal method that do not exceed the value β 0 √ max(ε r max ⊥ , ε r max // ) max(µ r max ⊥ , µ r max // ) for the phase constants (β), where β 0 , ε r max ⊥ (µ r max ⊥ ) and ε r max // (µ r max // ) indicate the respective phase constant in free space, the maximum transverse, and the maximum longitudinal relative permittivities (permeabilities) among the layers of the considered line.

FORMULATION OF THE PROPOSED FAST NUMERICAL INTEGRAL APPROACH INVOLVING NONUNIFORM DISCRETIZATION STEP SIZE
Solving the homogeneous system requires calculating the transverse inner products of basis functions with trial functions of the dispersion matrix, deducted from the integral formulation resulting from Eq. (6); the integration boundaries are given by the width of the slot (or metallic strip) related to the considered planar line. Nevertheless, because cosine trial functions with edge effects cannot be analytically integrable [17], the related inner products should be determined numerically. Then, the integral computation can be further improved by using a non-uniformly distributed discretization step with sinusoidal variations. Because of the shape of these trial functions, the integral discretization step needs to be large in the middle of the slot(s) (or of the metallic strip(s)) and very small close to the metallic strip edges, as shown in Fig. 4.  Here, x 0 and x N a indicate the respective first and last abscissa (or the first and the last position in the x-axis), involved in the integration boundaries of the inner products, while "i" and N a indicate the respective ith discretization step and the number of elementary areas of the numerical integration in the insulating (or metallic) region.
Furthermore, if the origin of x-axis is chosen on the top of the transmission planar lines such as in all asymmetrical structures as well as in symmetrical coplanar and coupled microstrip type structures, the ith discretization step can be expressed as: knowing that i = {0, 1, 2, . . . , N a} and x i ∈ [x 0 , x N a ]. However, by taking 1 ≤ i ≤ N a, the width of each elementary area is determined by: which represents also the nonuniform discretization parameter of the numerical integration. Moreover, according to Fig. 5, when the origin is chosen in the middle of the analyzed symmetrical finline or microstrip structure, x i and H i are expressed as follows: Remind that because of the symmetry of the proposed structures cited above, only half of the transverse section was considered, involving a reduction of the CPU time by 2. In the case of uniform discretization (i.e., all elementary areas have the same width H i = H), Eqs. (15), (17) and (16), (18)  can be reduced to To avoid heavy numerical computation of the inner products involving elementary areas, sinusoidal trial functions with no edge effects can be used, since they are analytically integrable. Such functions can be obtained by only considering the numerator of every cosine trial functions reported in Tables 1 and 2. Note that, for symmetrical and asymmetrical structures, such sinusoidal functions with no edge effects satisfy the boundary and proportionality conditions (as well as the parity condition for the symmetrical case). Also, although these functions neither generate spurious solutions [22] nor involve complex calculations, the convergence of the solutions (obtained by the zero determinant condition) is slower with regard to the size of the dispersion matrix. In fact, improving the non-trivial solution convergence of the previous homogeneous system requires using a dispersion matrix of relatively larger size.

ANALYTICAL APPROACH FOR COMPUTING INNER PRODUCTS USING TRIAL FUNCTIONS WITH METALLIC EDGE EFFECTS
In this section, we focus on the expressions required to numerically compute the transverse inner products involving cosine trial functions with edge effects (regarding the inner products involving sinusoidal trial functions, they can be determined analytically because they are more easily integrable). To numerically integrate them, let us use trapezoidal or quadratic elementary areas.

Uniform Width of Elementary Areas
In the case of uniform discretization step size, while considering the origin on the top of the modeled structure, the well-known expressions of the inner products are [27]: • Quadratic Areas

Nonuniform Width Of Elementary Areas
As for the integration involving nonuniform areas, original expressions of the inner products are introduced below. To the best of the authors' knowledge, such innovative expressions have never been proposed. In addition, to further increase the efficiency of the proposed approach, a technique to accelerate the convergence with regard to the number of areas will be investigated in the next section. When the origin is taken on the top of the considered structure, the inner products, which depend only on the x variable, can be described by the following new expressions involving a nonuniform discretization step: • Quadratic Areas It should be noted that if the origin is taken in the middle of a given structure such as symmetrical finlines (or microstrip lines), it suffices to multiply the four above formulas by two, where x 0 = 0 and N a indicates for the number of elementary areas defined in only half of the insulating (or metallic) region.

CONVERGENCE STUDY WITH REGARD TO THE NUMBER OF AREAS
This section highlights the convergence of the inner products. It involves only cosine trial functions, with regard to the number of areas while introducing a process to accelerate this convergence. To this aim, let us consider a symmetrical finline and its dual configuration, i.e., the symmetrical suspended microstrip line. As shown in Fig. 6, it includes two layers filled with air (µ r1 = µ r3 = 1; ε r1 = ε r3 = 1) and one layer constituting the substrate of ε r2 = ε r = 2.     The position of the last discretization of the numerical integration will be first varied. So, let ∆ be the narrowest width of the last nonuniform elementary area in the slot or on the metallic strip, which is the closest to the metallic edge, as shown in Fig. 7. Its expression is deducted from Eq. (18), as The latter discretization is located by the abscissa x N a = w/2 − q∆ with respect to the metal edge while varying the factor q from 0 to 0.4. In fact, this last discretization is shifted by q∆ with respect to the metallic edge. In this case, the integration boundaries of these inner products are taken between 0 and w/2 − q∆.  Figure 8 illustrates the evolution of the inner products versus the number of nonuniform trapezoidal and quadratic areas for different positions of the last discretization (depending on the factor q), with w/a = 0.9, m = 2, and n = 100 for the symmetrical finline and n = 99 for the symmetrical suspended microstrip line.
The above figure shows a fast convergence of inner products versus the number of trapezoidal areas when the last discretization is shifted by 0.25∆ from the metallic edge (i.e., x N a = w/2 − 0.25∆). Therefore, this position of "x N a " can be considered to be adequate for fast and accurate analysis of planar structures. On the other hand, varying the position of the last discretization (by varying "q") has no influence on the convergence when quadratic areas are used.
To further emphasize the optimal position of the last discretization with respect to the convergence speed of the inner products, we report in Tables 3 and 4 the convergence of the inner products involving cosine trial functions as function of the number of elementary areas (Note that the trapezoidal form is involved only for finline and microstrip structures). Table 3. Convergence of inner products versus trapezoidal areas for finline structure.

Number of areas
Last position  These two tables clearly show that the convergence of the inner products becomes faster when the last discretization step is shifted by 0.25∆ from the metallic edge, knowing that x N a = w/2 − 0.25∆ represents the upper boundary of the numerical integration of the inner products. Nevertheless, to further verify the efficiency of the proposed nonuniform numerical integral technique compared to the uniform one, we displayed the evolution of the ox-and oz-inner products of the symmetrical finline ( Fig. 9(a)) and suspended microstrip line ( Fig. 9(b)), respectively, as function of the number of trapezoidal or quadratic areas for these two types of discretization steps. In this part, we considered two values for x N a namely, when the last discretization line is on the metallic edge (x N a = w/2) and when it is shifted by 0.25∆ from the latter (x N a = w/2 − 0.25∆), with w/a = 0.9, K = 2, N = 100 for the finline and N = 99 for the suspended microstrip line.
From this figure, we can deduce that when the last discretization is on the edge (x N a = w/2), the convergence was very slow requiring at least 1,000,000 and 3,000 elementary areas while using uniform and nonuniform discretization steps, respectively. However, and as expected, the convergence was considerably accelerated when the last discretization was shifted by 0.25∆ from the metallic edge (x N a = w/2 − 0.25∆) implying more reduction of trapezoidal or quadratic areas. Figure 10 shows a comparative study of the determinant of the dispersion matrix obtained from the two numerical integration methods using uniform/nonuniform discretization steps, when trial functions with edge effects are considered. We note a faster convergence according to the number of nonuniform trapezoidal areas while shifting the last discretization by 0.25∆ from the metallic edge, knowing that this convergence does not exceed 500 areas.
According to the obtained results, the computational integral method using nonuniform trapezoidal areas can be an appropriate choice for proper and reliable analysis of planar transmission lines. Therefore, we report in Fig. 11 the convergence of the phase constant (β) of the dominant mode with regard to the nonuniform trapezoidal areas for the two above positions of the last discretization (x N a ) and by taking at least 2000 modes and 2 trial functions per component. Then, for a number not exceeding 400 nonuniform trapezoidal areas, we reached the convergence for x N a = w/2 − 0.25∆. So, we considered only nonuniform trapezoidal areas when the last discretization was located at w/2 − 0.25∆. Fig. 12 shows the convergence of the dominant mode phase constant (β) with respect to the number of nonuniform trapezoidal areas, for different values of normalized widths w/a of the slot (metallic strip) of the symmetrical finline (suspended microstrip line). Figure 12 shows a slower convergence when the slot (metallic strip) for the symmetrical finline (suspended microstrip line) is wider. Note that, because we used a sinusoidal nonuniform discretization step, this convergence does not exceed 500 areas whatever the case. Figure 13 highlights the convergence with regard to the number of trial functions (K) and modes (N ) while choosing narrow or wide slot (or metallic strip). In this case, both trial functions with and without edge effects were used. We can notice that using trial functions with edge effects results in both a slower convergence with regard to N and a faster one with regard to K, requiring up to 3000 modes and between 2 and 8 trial functions. On the other hand, by using trial function without edge effects, the convergence is slower with respect to K, requiring at least 20 trial functions per component to reach a stable solution, thus implying the use of a relatively large size dispersion matrix, knowing that this convergence is obtained from 300 modes.
To demonstrate the reliability of the proposed numerical integral technique, we evaluated β versus the normalized width while selecting either 4 trial functions by component when edge effects are considered or 20 when edge effects are disregarded (Fig. 14). We accordingly noted a close similarity between the results that we obtained from these two types of trial functions, thus validating the numerical integral approach based on a proper computation of the inner products. In fact, the use of our numerical approach with cosine trial functions taking into account edge effects leads to the same results derived from an analytical integral calculation based on cosine trial functions without edge effects.
Nevertheless, the same figure shows a decrease of β versus the normalized slot width of the finlinetype structure, and a decrease then an increase of β versus the normalized metallic strip width of the microstrip-type structure. In fact, this increase is due to the influence of the shielding vertical walls on the electromagnetic field distribution, especially when the edges of the metallic strip are closer to these walls.
The efficiency of the proposed modeling method depends not only on the accuracy of the expected results but also on the required computing time. In fact, CPU time is another key criterion to evaluate the performance of our approach. As shown in Fig. 15 Indeed, the proposed approach helps reducing the CPU time by calculating the determinant of the dispersion matrix for only one value of the phase constant. Also, this figure clearly shows the influence of the dispersion matrix size, which ascertains a suitable compromise between memory size, accuracy, and CPU time. From this figure, we can also deduce that the numerical integration approach involving trapezoidal areas reduces the CPU time by a factor of 1.5 compared to the case where quadratic areas are used, since the mathematical expression of the numerical integration using trapezoidal areas is simpler than the one involving quadratic elementary areas. As illustration, for K = 8 and N = 3000, the CPU time is about 20 s and 30 s for respectively trapezoidal and quadratic areas. So, the trapezoidal approach is about 1.5 times faster than the quadratic approach. We, therefore, retained the trapezoidal numerical approach including a nonuniform mesh when trial functions with edge effects are selected (while shifting the last discretization step by 0.25∆ from the metallic edge).

NUMERICAL EXAMPLES INCLUDING SYMMETRICAL/ASYMMETRICAL STRUCTURES
To further confirm the efficiency of the proposed numerical integral approach, the obtained results have been compared to those published in the technical literature, by using a sufficient number of areas, modes, and trial functions with and without edge effects. For this purpose, we derived the effective permittivity, normalized wavelength, and phase constant for the dominant mode as function of different physical and electrical parameters of the retained symmetrical/asymmetrical structures (where the study region is taken either in the slot(s) or on the metallic strip(s), according to the type of considered lines).

Finline-Type Structure
Starting with symmetrical finline structures with isotropic/anisotropic substrates, we plotted in Fig. 16 the effective permittivity vs. the normalized slot (w/a) where the two types of trial functions (with and without edge effects) were used. The retained anisotropic materials, for such finline, are Epsilam 10, sapphire, and boron nitride characterized by their electrical tensorsε r = [13, 10.2, 13],ε r = [9.4, 11.6, 9.4] andε r = [5.12, 3.4, 5.12], respectively. Note that the decrease of the effective permittivity with respect to the w/a ratio is in close agreement with the results published by Shalaby and Kumar [28] with an overall average relative error less than 1.2%. Note that the two types of trial functions yield practically the same results by selecting a sufficient number of modes (N = 3000), areas for trial functions with edge effects (N a ≤ 700), and trial functions by component (i.e., K ≤ 8 and K ≤ 25 for trial functions with and without edge effects, respectively).  Figure 17 plots the frequency variation of the dispersion parameters of the microstrip-type structure with a maximum of 3 layers such as symmetrical microstrip with overlay. Note that the structure consists of a metallic strip inserted between two identical isotropic (anisotropic) dielectric layers of

Microstrip-Type Structure
ε r = 4.08 (ε r = [5.12-3.4-5.12]). We can observe an increase of the effective permittivity (ε eff ) versus both the w/λ 0 ratio and the height of the 2nd layer (b 2 ) of the covered symmetrical microstrip, with λ 0 being the free space wavelength.
On the other hand, for b 2 /w = 1, the frequency has nearly no effect on dispersion since ε eff is almost constant, because this structure becomes a homogeneous stripline (b 1 = 0), consisting of a single medium in which the same concentration of the EM fields is confined. So, by using trial functions with and without edge effects, our results generally agree well with those of Marques and Horno [29]. Nevertheless, when the strip is very narrow with respect to the wavelength (w/λ 0 < 0.02), trial functions without metallic edge effects led to less accurate values than those obtained from trial functions that include edge effects. Hence, in this case, favoring the use of trial functions with edge effect ensures better accuracy of the dispersion parameters. Moreover, by referring to the work of Marques and Horno who used the spectral method including basic functions with edge effects [29], we obtained close results with an average error of 0.32% and 0.44% while using trial functions with and without edge effects, respectively.

Coupled Microstrip-Type Structure
The next step is to consider coupled microstrip-type structures. Fig. 18 illustrates the quasi-static analysis by varying the effective permittivity versus physical parameters of symmetrical/asymmetrical coupled microstrip-type structures. In Fig. 18(a), we notice an increase of the even-and odd-mode effective permittivity (ε eff ) versus the b 3 /b 2 ratio for the shielded inverted coupled microstrip lines with isotropic/anisotropic substrate. In fact, by increasing the height b 3 , the EM fields are less attracted towards the lower horizontal wall of the shielding, thus increasing the concentration of these fields in the substrate and, therefore, the effective permittivity. However, when the metallic strips are sufficiently away from the lower horizontal wall by further increasing the height b 3 , this latter wall will have no effect on the propagation characteristics, since most of the EM fields will be attracted into the substrate, thus stabilizing ε eff when the ratio b 3 /b 2 increases. Furthermore, this figure highlights a close agreement with results published in [30] with an average relative error that does not exceed 0.84% when the two types of trial functions were used. width (d 2 /d 1 ) of the metallic strips, this figure highlights a decrease of ε eff for both suspended coupled lines and coupled lines with an overlay while it is increasing for coupled microstrip lines. Moreover, according to this figure, when d 2 /d 1 < 1, trial functions without edge effects lead to less accurate values, which favors the use of trial functions with edge effects. Indeed, the total average relative errors resulting from trial functions with and without edge effects are 0.16% and 0.63%, respectively, while comparing our values with those of Kitazawa and Mittra [31].

Coplanar-Type Structure
After selecting a sufficient number of modes, trapezoidal areas, and trial functions, we computed the normalized phase constant and effective permittivity of shielded coplanar structures with isotropic magnetic/non-magnetic substrates and inserted in a WR28-type waveguide of section 3.556×7.112 mm 2 .
To this aim, we first considered a symmetrical coplanar waveguide printed on a magnetic isotropic substrate while varying the values of the relative permeability. In this case, as shown in Fig. 19(a), there is a clear influence of the magnetic isotropy on the dispersion curve by plotting the square of the normalized phase constant (β/β 0 ) 2 of the dominant mode (even mode) versus frequency. Note that this ratio was numerically determined by Mazé et al. [32] who defined it as the product of the effective permittivity with the effective permeability (ε eff * µ eff ). In addition, the dispersion becomes higher with the increase of µ r , with a significant influence of the non-magnetic isotropy on the dispersion characteristics. On the other hand, when µ r = 1 (non-magnetic isotropy case), a very slight dispersion is observed. With an average relative error not exceeding 0.92%, our computed results can be successfully compared to those reported in [32], knowing that the two types of selected trial functions led to similar results. Figure 19(b) shows the evolution of the effective permittivity of the shielded asymmetrical coplanar structure versus the width of the lateral metallic strip (C 1 ) for both C-and π-modes. From this figure, we see a quasi-stability of ε eff with respect to C 1 for both C-and π-modes, highlighting a negligible influence of the vertical metallic walls on the dispersion characteristics when these walls are sufficiently away from the central metallic strip. For both trial functions with and without edge effects, the obtained results are in good agreement with those reported by Schmidt et al. [33] with an average error less than 0.35%.

CONCLUSION
An integral method based on the EM mathematical concept of operators was proposed to analyze arbitrary isotropic/anisotropic multilayered structures, involving the Galerkin's procedure in the modal domain. For this purpose, two types of trial functions were selected namely, sinusoidal trial functions, which allow lighter analytical calculations but with slower convergence (with respect to the dispersion matrix size), and trial functions with metallic edge effects that favor a faster convergence by using numerical integration with a nonuniform discretization step, but at the expense of higher mathematical development. Furthermore, optimal convergence of the propagation parameters was investigated leading to close agreement of the obtained results with those published in the technical literature. We should notice that in some cases, i.e., when the width of the slot or metallic strip is much smaller than that of the concerned structure, trial functions with edge effects will improve the accuracy. Thus, this numerical integral approach constitutes an advanced modeling tool for the efficient and reliable design of RF/microwave planar circuits.