Partial-inductance retarded partial coefficients: Their exact computation based on the Cagniard-DeHoop technique

Engineering Analysis with Boundary Elements 149 (2023) 86–91 Available online 20 January 2023 0955-7997/© 2023 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound Partial-inductance retarded partial coefficients: Their exact computation based on the Cagniard–DeHoop technique Martin Stumpfa, Fabrizio Loretob,∗, Giuseppe Pettanicec, Giulio Antoninib aLerch Laboratory of EM Research, Department Radio Electronics, FEEC, Brno University of Technology, Brno, Czech Republic bDepartment of Industrial and Information Engineering and Economics, Università degli Studi dell’Aquila, L’Aquila, Italy cDepartment of Engineering and Information Sciences and Mathematics, Università degli Studi dell’Aquila, L’Aquila, Italy ARTICLE INFO Keywords: Partial Element Equivalent Circuit method Cagniard–DeHoop technique Computational electromagnetics Partial inductances Time-domain modeling ABSTRACT The Partial Element Equivalent Circuit (PEEC) method is a well recognized integral-equation (IE) technique to solve Maxwell’s equations. Similarly to the method of moments (MoM), the electromagnetic (EM) interactions between currents and between charges are described in terms of integrals. In contrast to the standard MoM, the PEEC method keeps the electric and magnetic coupling phenomena separate, which leads to different interaction integrals to be computed. These integrals admit simplified solutions for the case of the static free-space Green’s function and orthogonal geometries but their applicability is limited to electrically small problems only. When the full-wave free-space Green’s function is considered, the integrals are typically computed in the frequency domain (FD) by resorting to Gaussian quadrature schemes. The accuracy and efficiency of such schemes is a delicate issue. Therefore, recent works have investigated the possibility of applying the Cagniard–DeHoop (CdH) technique to calculate the interaction integrals for zero-thickness elementary domains. In this paper, we close the loop and shall apply the CdH technique to calculate the partial-inductance between two elementary bricks as prescribed by the PEEC technique exactly in the time domain (TD). The analytical approach is demonstrated on the interaction between two bricks as it occurs in the modeling of the magnetic field coupling between volumetric currents. The accuracy of the proposed approach is (successfully) tested for two representative cases. 1. Introduction The notion of partial inductance, as introduced by Dr. A. E. Ruehli in 1972 [1], is a fundamental concept on which the PEEC method lays its foundations. The PEEC method is an IE method capable of analyzing an EM scattering problem by means of an equivalent circuit representation [2]. The PEEC method can be formulated with the aid of standard EM contrast-source integral representations [3, Sec. 28.9], in which the Green’s functions apply to the (typically homogeneous, isotropic and loss-free) medium of the embedding. Hence, the equivalent circuit can be directly associated with a discretized scatterer, where elementary volumes and surfaces are assumed to radiate in the background medium (e.g. the free-space). The PEEC method has been primarily employed in the solution of FD problems and, in this context, interactions integrals have been computed either making quasi-static assumptions (e.g., neglecting the retardation terms in the integrals) [4] or by resorting to quadrature schemes in the frequency domain (FD) [5] or by means of the Taylor expansion of the Green’s function [6]. Similar lines of reasoning are commonly followed to extract the singularity ∗Corresponding author. E-mail address: [email protected] (F. Loreto). of Green’s function in standard integral-equation formulations (e.g., [7– 9]). Owing to the ever increasing interest in accurate TD simulations in the field of Electromagnetic Compatibility (EMC), however, fullwave TD formulations of the partial inductance are becoming more and more important. Accordingly, this issue has been recently addressed analytically for relatively simple 2-D zero-thickness patches [10,11]. In this article, we shall apply the CdH technique [12] to calculate a partial-inductance retarded PEEC coefficient exactly in the TD. The presented analytical results apply to the interaction between two bricks (= right parallelepipeds). The CdH technique is a joint-transform method that has been originally developed to analytically analyze the seismic-wave propagation in horizontally layered media (e.g. [13, 14]). More recently, it has been demonstrated that this sophisticated inversion methodology can also be useful for constructing purely numerical solutions. Indeed, the CdH technique is a key ingredient in the Cagniard–DeHoop method of moments, a novel TD-IE technique for the TD analysis of EM radiation and scattering problems [15,16]. Moreover, our initial studies (see [10,11]) analyzing zero-thickness PEEC https://doi.org/10.1016/j.enganabound.2023.01.008 Received 9 August 2022; Received in revised form 8 January 2023; Accepted 8 January 2023 Engineering Analysis with Boundary Elements 149 (2023) 86–91 87 M. Stumpf et al. Fig. 1. Two interacting brick elements. coefficients have demonstrated that the CdH inversion is a promising strategy for achieving their analytical expressions in the TD. Accordingly, this paper reports on the recent results of our efforts to develop a quasi-static-approximation-free PEEC solver that is equipped with exact, CdH-based TD volumetric partial-inductance. It is anticipated that the presented approach can be employed along with the precorrected fast Fourier transform (FFT) [17] and the FFT-based approach to accelerate the matrix–vector products [18–20]. The paper is organized as follows: in Section 2the problem formulation is presented. Consequently, the TD analytical solution is given in Section 3, which is supplemented with Appendix. Finally, in Section 4, three numerical applications are presented and successfully validated. 2. Problem formulation A PEEC model is represented through a set of partial elements, the value of which is found upon evaluating spatial integrals over the surfaces/volumes of interacting discretization elements. In this work we shall analyze the interaction of two bricks (= right parallelepipeds) (see Fig. 1). In particular, with reference to [21, Eq. (4)], we shall study a retarded partial potential coefficient expressed through a double integral  𝐿𝑚𝑛(𝑠) = 𝜇0 𝑚𝑛∫𝒓∈𝑚 d𝑉∫𝒓′∈𝑛 𝑔(𝒓−𝒓′, 𝑠)d𝑉′,(1) where 𝑚= {−𝛥𝑚 𝑥∕2 < 𝑥 −𝑥𝑚< 𝛥𝑚 𝑥∕2,−𝛥𝑚 𝑦∕2 < 𝑦 −𝑦𝑚< 𝛥𝑚 𝑦∕2,−𝛥𝑚 𝑧∕2 < 𝑧−𝑧𝑚< 𝛥𝑚 𝑧∕2} and 𝑛= {−𝛥𝑛 𝑥∕2 < 𝑥 −𝑥𝑛< 𝛥𝑛 𝑥∕2,−𝛥𝑛 𝑦∕2 < 𝑦 −𝑦𝑛< 𝛥𝑛 𝑦∕2,−𝛥𝑛 𝑧∕2 < 𝑧 −𝑧𝑛< 𝛥𝑛 𝑧∕2}, where 𝛥𝑚,𝑛 𝑥>0,𝛥𝑚,𝑛 𝑦>0and 𝛥𝑚,𝑛 𝑧>0 denote the spatial discretization steps in the 𝑥-, 𝑦and 𝑧-direction, respectively. Furthermore, 𝑠is the Laplace-transform parameter with Re(𝑠)>0, and 𝑚,𝑛 are cross sections of the volumes 𝑚,𝑛, respectively, that are perpendicular to the corresponding electric-current flows. Next, 𝑔(𝒓, 𝑠) = exp(−𝑠|𝒓|∕𝑐) 4𝜋|𝒓|(2) is the free-space Green’s function of the 3-D scalar modified Helmholtz equation and 𝑐= (𝜀𝜇)−1∕2 >0denotes the pertinent (real-valued and positive) EM wave speed. 3. Problem solution The retarded partial potential coefficient as expressed through Eq. (1) will be next transformed to the TD analytically with the aid of the CdH technique. Pursuing this approach and assuming that Fig. 2. Two interacting cube elements. |𝑧𝑚−𝑧𝑛|>(𝛥𝑚 𝑧+𝛥𝑛 𝑧)∕2, one may express the TD original of Eq. (1), further denoted by 𝐿𝑚𝑛(𝑡)(in henry/second = ohm), as follows: 𝐿𝑚𝑛(𝑡)=(𝜇0∕𝑚𝑛)[𝐽(|𝑧𝑚−𝑧𝑛|+𝛥𝑚𝑛+ 𝑧, 𝑡) −𝐽(|𝑧𝑚−𝑧𝑛|+𝛥𝑚𝑛− 𝑧, 𝑡) − 𝐽(|𝑧𝑚−𝑧𝑛|−𝛥𝑚𝑛− 𝑧, 𝑡) +𝐽(|𝑧𝑚−𝑧𝑛|−𝛥𝑚𝑛+ 𝑧, 𝑡)],(3) where 𝐽(𝑧, 𝑡) = 𝐼(𝑥𝑚−𝑥𝑛+𝛥𝑚𝑛+ 𝑥, 𝑧, 𝑡) −𝐼(𝑥𝑚−𝑥𝑛+𝛥𝑚𝑛− 𝑥, 𝑧, 𝑡) − 𝐼(𝑥𝑚−𝑥𝑛−𝛥𝑚𝑛− 𝑥, 𝑧, 𝑡) +𝐼(𝑥𝑚−𝑥𝑛−𝛥𝑚𝑛+ 𝑥, 𝑧, 𝑡),(4) and 𝐼(𝑥, 𝑧, 𝑡) = 𝐾(𝑥, 𝑦𝑚−𝑦𝑛+𝛥𝑚𝑛+ 𝑦, 𝑧, 𝑡) −𝐾(𝑥, 𝑦𝑚−𝑦𝑛+𝛥𝑚𝑛− 𝑦, 𝑧, 𝑡) −𝐾(𝑥, 𝑦𝑚−𝑦𝑛−𝛥𝑚𝑛− 𝑦, 𝑧, 𝑡) +𝐾(𝑥, 𝑦𝑚−𝑦𝑛−𝛥𝑚𝑛+ 𝑦, 𝑧, 𝑡),(5) where we used 𝛥𝑚𝑛± 𝑥,𝑦,𝑧 = (𝛥𝑚 𝑥,𝑦,𝑧 ±𝛥𝑛 𝑥,𝑦,𝑧)∕2,(6) respectively. Here, 𝐾(𝑥, 𝑦, 𝑧, 𝑡)represents the TD original of the generic slowness integral, the definition and inversion of which is presented in Appendix. Finally, we emphasize that Eq. (3) applies to the configuration where |𝑧𝑚−𝑧𝑛|> 𝛥𝑚𝑛+ 𝑧. The case |𝑧𝑚−𝑧𝑛|< 𝛥𝑚𝑛+ 𝑧must be analyzed separately. 4. Numerical examples In this section, numerical examples related to three different geometries are presented. For validation purposes, the results obtained through the proposed method are compared with those obtained through the numerical-inversion of the Laplace transform (NILT) approach [22–24]. 4.1. Two interacting cubes As a particular application, the resulting TD expression (3) has been implemented in MATLAB®and applied to the case of two interacting cubes (see Fig. 2). The evaluations are performed in the finite time window {0 ≤𝑐𝑡∕𝑅𝑚𝑛 ≤2}, where 𝑅𝑚𝑛 = [(𝑥𝑚−𝑥𝑛)2+ (𝑦𝑚−𝑦𝑛)2+ (𝑧𝑚− 𝑧𝑛)2]1∕2 represents the center-to-center distance between two (identical) cubes located at •(𝑥𝑚, 𝑦𝑚, 𝑧𝑚) = (0,0,0), Engineering Analysis with Boundary Elements 149 (2023) 86–91 88 M. Stumpf et al. Fig. 3. The TD coefficient for the two identical cubes: comparison between the proposed technique and the dNILT2 - Hermite technique.. Fig. 4. Spectrum comparison for the interaction between two identical cubes. •(𝑥𝑛, 𝑦𝑛, 𝑧𝑛) = (2𝛥𝑥,2𝛥𝑥,2𝛥𝑥). In the present example we take 𝛥𝑥=𝛥𝑚,𝑛 𝑥=𝛥𝑚,𝑛 𝑦=𝛥𝑚,𝑛 𝑧= 1.0 mm. The resulting pulse shape of the TD coefficient is shown in 3. Here, for the sake of validation, the results obtained through the proposed CdHbased and the (referential) dNILT2 - Hermite technique are presented. For a detailed description of the referential methodology we refer the reader to [24]. The circle points are the initial samples needed for the interpolant building. Finally, the FD counterpart of the two TD responses is sketched in Fig. 4. As can be seen, the computed results show good correspondence with dissimilarities occurring from 400 GHz up. These discrepancies are, however, virtually negligible in technical applications. 4.2. Interactions of a system of cubes As a further example we consider the geometry depicted in Fig. 5, where the mutual partial inductance between the cube with the center located at the axes origin and the others are considered. All the cubes have sides 𝛥𝑥=𝛥𝑦=𝛥𝑧= 1.0 mm. The cubes system is composed by 27 elements, three for each dimension. The center coordinates span a range [2𝛥𝑥− 8𝛥𝑥]. In Fig. 6 are sketched all the TD coefficients related to the system, each one computed through the proposed technique and compared to the dNILT2 - Hermite technique [24]. The latter Fig. 5. Cubes system geometry. Fig. 6. TD coefficients for the cubes system geometry. technique is based on the construction of an accurate interpolator, starting from the knowledge of the values of the function and its first higher order derivatives at the starting points, encircled in red in the figure (generally up to the fifth order is sufficient). 4.3. Induced voltage on a cube by the currents flowing in four bricks In Fig. 7 are shown four identical parallelepipeds: 1,2,3,4, with sides: 𝛥𝑥= 1.5 mm, 𝛥𝑦= 0.5 mm, 𝛥𝑧= 0.25 mm, above a cube 0, with sides: 𝛥0 𝑥=𝛥0 𝑦=𝛥0 𝑧= 1.0 mm, and center located at the axes origin. The parallelepipeds are considered at the same height and are collocated unsymmetrically in the 𝑥–𝑦plane, with respect to the cube. In particular, the coordinates of the parallelepipeds are: •(𝑥1, 𝑦1, 𝑧1) = (−2.5 mm,0.5 mm,3𝛥0 𝑥), •(𝑥2, 𝑦2, 𝑧2) = (−2.5 mm,5 mm,3𝛥0 𝑥), •(𝑥3, 𝑦3, 𝑧3) = (11 mm,0.5 mm,3𝛥0 𝑥), •(𝑥4, 𝑦4, 𝑧4) = (11 mm,5 mm,3𝛥0 𝑥). The overall induced voltage on the cube 0by the system of parallelepipeds can be computed through a combination of four convolution integrals as: 𝑣𝐿0(𝑡) = 4 ∑ 𝑛=1 ∫𝑡 0 𝐿𝑝0,𝑛 (𝑡−𝜏)d𝑖𝑛(𝜏) d𝜏d𝜏(7) Engineering Analysis with Boundary Elements 149 (2023) 86–91 89 M. Stumpf et al. Fig. 7. Geometry for the computation of the induced voltage on a cube by the currents flowing in the bricks system. Fig. 8. Induced voltage on the cube by the currents flowing in the system of bricks. where 𝑖𝑛(𝑡)is the impressed current on each parallelepiped, 𝑛= 1, ⋯,4, flowing in the 𝑥direction. The impressed current 𝑖𝑛(𝑡), for each parallelepiped, is assumed to exhibit a windowed-power (WP) waveform [25]: 𝑖𝑛(𝜏, 𝑡) = 𝑡′𝜏(2 − 𝑡′)𝜏H(𝑡′)H(2 − 𝑡′)(8) where H(𝑡)is the Heaviside unit-step function (H(𝑡)=0if 𝑡 < 0, H(0) = 1∕2,H(𝑡) = 1 if 𝑡 > 0), 𝑡′=𝑡∕𝑡r,𝑡rbeing the pulse rise time. We choose 𝜏= 2 and 𝑡r= 4.6ps. The induced voltage on the cube 0is depicted in Fig. 8, where it is observed an excellent agreement between the CdH method and the NILT-based method. 5. Conclusions The PEEC method requires that interaction integrals describing the magnetic field coupling between elementary volumetric regions be computed, namely partial inductances. In the frequency domain, this is usually done by resorting to quadrature schemes. In the TD, the use of over-simplifying assumptions on the propagation delay leads to approximate results with a negative impact on physical properties such as the causality and stability of the model. In this work, quasiclosed-form for TD retarded partial inductances have been derived using the Cagniard–DeHoop (CdH) technique. A pertinent integration path deformation in the complex slowness plane allows to obtain semianalytical forms of the transient interaction integrals for a pair of orthogonal bricks as they occur in the PEEC method using Manhattantype meshes or voxellization techniques. The proposed approach has been tested for representative test cases by comparison with other numerical methods, always exhibiting a very good agreement. Fig. 9. Complex slowness planes. (a) 𝜎-plane with the CdH-path for 𝑦 < 0; (b) 𝜅-plane with the CdH-paths for 𝑥 < 0. Declaration of competing interest The authors declare that they have no conflict of interest. Data availability The presented data are available upon request from the authors. Acknowledgments The research of Martin Stumpf was supported by the Czech Science Foundation under Grant No. 20-01090S. Appendix. The generic integral The integral representation to be transformed to TD has the following form  𝐾(𝑥, 𝑦, 𝑧, 𝑠) = (𝑠 2i𝜋)2∫𝜅∈K0 exp(𝑠𝜅𝑥) 𝑠2𝜅2d𝜅 ×∫𝜎∈S0 exp{−𝑠[−𝜎𝑦 +𝛤(𝜅, 𝜎)𝑧]} 𝑠2𝜎2 d𝜎 2𝑠3𝛤3(𝜅, 𝜎)(9) for 𝑥∈R,𝑦∈R,{𝑧∈R;𝑧≥0} and {𝑠∈R;𝑠 > 0}, where K0and S0are the integration paths extending along Re(𝜅)=0and Re(𝜎)=0, respectively, that are indented to the right with semi-circular arcs with centers at the origins and vanishingly small radii (see Fig. 9). Finally, 𝛤(𝜅, 𝜎), being the slowness parameter along the 𝑧-direction, is defined as 𝛤(𝜅, 𝜎) = (1∕𝑐2−𝜅2−𝜎2)1∕2 with Re(𝛤)≥0.(10) The generic integral will next be transformed to the TD with the aid of the CdH technique (see [12] and [16, Ch. 2]). To that end, the integration contour in the complex 𝜎-plane, S0, is by virtue of Jordan’s lemma and Cauchy’s theorem [3, p. 1054] deformed into a CdH path, say ∪∗(here ∗denotes the complex conjugate), along Engineering Analysis with Boundary Elements 149 (2023) 86–91 90 M. Stumpf et al. which −𝜎𝑦 +𝛤(𝜅, 𝜎)𝑧=𝑢𝑑𝛺(𝜅)for {1 ≤𝑢 < ∞} with 𝑑2=𝑦2+𝑧2 and 𝛺(𝜅) = (1∕𝑐2−𝜅2)1∕2 is satisfied (see Fig. 9a). Upon combining the contributions from and ∗, the inner integral with respect to 𝜎can be cast into the integral with respect to the (real-valued and positive) parameters 𝑢. In addition, the contribution from the (double) pole singularity at 𝜎= 0 must be for 𝑦 > 0accounted for. The thus expressed inner integral is subsequently substituted back in Eq. (9), which yields  𝐾(𝑥, 𝑦, 𝑧, 𝑠) =  𝑀(𝑥, 𝑦, 𝑧, 𝑠) +  𝑁(𝑥, 𝑦, 𝑧, 𝑠),(11) where  𝑀=1 2𝜋i 𝑑4 2𝜋𝑠3∫∞ 𝑢=1 𝑦2𝑧2−𝑢2(𝑢2− 1)(𝑦4− 6𝑦2𝑧2+𝑧4) (𝑢2𝑑2−𝑦2)2(𝑢2𝑑2−𝑧2)2 ×d𝑢 (𝑢2− 1)1∕2 ∫𝜅∈K0 exp{−𝑠[−𝜅𝑥 +𝛺(𝜅)𝑢𝑑]} d𝜅 𝑠2𝜅2𝛺4(𝜅)(12) and  𝑁=1 2𝜋i 𝑦H(𝑦) 2𝑠2∫𝜅∈K0 exp{−𝑠[−𝜅𝑥 +𝛺(𝜅)𝑧]} d𝜅 𝑠2𝜅2𝛺3(𝜅),(13) where H(𝑦)has again the meaning of the Heaviside unit-step function, i.e. H(𝑦)=0if 𝑦 < 0,H(0) = 1∕2,H(𝑦)=1if 𝑦 > 0. First, we shall describe the transformation of  𝑀as given by Eq. (12). For this purpose, the integration contour in the complex 𝜅-plane, K0, is deformed into a CdH path, say ∪∗, along which −𝜅𝑥+𝛺(𝜅)𝑢𝑑 =𝜏for {𝑅(𝑢)∕𝑐≤𝜏 < ∞} with 𝑅(𝑢)=(𝑥2+𝑢2𝑑2)1∕2 >0is satisfied (see Fig. 9b). In the resulting expression, we combine the contributions from and ∗and change the order of the integrations according to, symbolically ∫∞ 𝑢=1 d𝑢∫∞ 𝜏=𝑅(𝑢)∕𝑐 d𝜏→∫∞ 𝜏=𝑅(1)∕𝑐 d𝜏∫𝑈(𝑐𝜏) 𝑢=1 d𝑢(14) where 𝑈(𝑐𝜏) = (𝑐2𝜏2∕𝑑2−𝑥2∕𝑑2)1∕2. Upon carrying out the integration with respect to 𝑢, Eq. (12) can be cast into the following form  𝑀=𝑐6 2𝜋2𝑠5∫∞ 𝜏=𝑅(1)∕𝑐 exp(−𝑠𝜏)(𝑥, 𝑦, 𝑧, 𝑐𝜏)d𝜏 + 𝑃(𝑥, 𝑦, 𝑧, 𝑠)(15) where  𝑃arises from the (double) pole singularity at 𝜅= 0. Both terms on the right-hand side of Eq. (15) have the form that allows their straightforward transform to the original domain. In this procedure, Lerch’s uniqueness theorem applying to the real-valued and positive Laplace-transform parameter is an essential result that we rely on [26, Appendix]. The transformation of  𝑁(see Eq. (13)) follows similar lines of reasoning. Indeed, the original integration contour, K0, is first replaced with a new CdH path along which −𝜅𝑥+𝛺(𝜅)𝑧=𝜏is met for all {𝜌∕𝑐≤ 𝜏 < ∞}, where 𝜌2=𝑥2+𝑧2. Combining again the contributions from the hyperbolic arcs in the lower and upper halves of the complex 𝜅-plane, we end up with an integral with respect to 𝜏that can be expressed as  𝑃(𝑦, 𝑥, 𝑧, 𝑠)(cf. Eq. (15)). Representing further the contribution from the (double) pole singularity at 𝜅= 0 by  𝑄(𝑥, 𝑦, 𝑧, 𝑠), we arrive at  𝑁= 𝑃(𝑦, 𝑥, 𝑧, 𝑠) +  𝑄(𝑥, 𝑦, 𝑧, 𝑠).(16) Upon substituting Eqs. (15) with (16) in (11) and transform the result to the TD, we finally get 𝐾(𝑥, 𝑦, 𝑧, 𝑡) = 𝑐 48𝜋2∫𝑐𝑡 𝑣=𝑅 (𝑐𝑡 −𝑣)4(𝑥, 𝑦, 𝑧, 𝑣)d𝑣 +𝑃(𝑥, 𝑦, 𝑧, 𝑡) + 𝑃(𝑦, 𝑥, 𝑧, 𝑡) + 𝑄(𝑥, 𝑦, 𝑧, 𝑡).(17) The function behind the integral sign is given by (𝑥, 𝑦, 𝑧, 𝑣) = ∫𝜋∕2 𝜓=0 𝑓(𝑥, 𝑦, 𝑧, 𝑣, 𝜓) ×𝑝2−𝑞2(𝑈2− 1) sin2(𝜓)[cos2(𝜓) + 𝑈2sin2(𝜓)] {(𝑧2∕𝑑2) cos2(𝜓) + [(𝑣2−𝑟2)∕𝑑2] sin2(𝜓)}2 ×d𝜓 {(𝑦2∕𝑑2) cos2(𝜓) + [(𝑣2−𝜌2)∕𝑑2] sin2(𝜓)}2,(18) with 𝑟2=𝑥2+𝑦2,𝑝2=𝑦2𝑧2∕𝑑4,𝑞2=𝑦4∕𝑑4− 6𝑝2+𝑧4∕𝑑4,𝑈2= 𝑣2∕𝑑2−𝑥2∕𝑑2and 𝑓=𝑣∕𝑑2 𝑈3[(𝑥2∕𝑑2) sin2(𝜓)+(𝑣2∕𝑑2− 1) cos2(𝜓)]2 ×{(3𝑥2 𝑑2+𝑣2 𝑑2)[cos2(𝜓) + 𝑈2sin2(𝜓)]3 +(4𝑥4 𝑑4−11𝑥2𝑣2 𝑑4−𝑣4 𝑑4)[cos2(𝜓) + 𝑈2sin2(𝜓)]2 −(𝑥6 𝑑6+5𝑥4𝑣2 𝑑6− 10 𝑥2𝑣4 𝑑6)[cos2(𝜓) + 𝑈2sin2(𝜓)] − 2 𝑥8 𝑑8+ 7 𝑥6𝑣2 𝑑8− 5 𝑥4𝑣4 𝑑8}.(19) The remaining terms in the final TD result (17) can be expressed as follows 𝑃(𝑥, 𝑦, 𝑧, 𝑡) = 𝑐𝑥𝑑3H(𝑥) 12𝜋{6|𝑦|𝑧𝑐𝑡 𝑑3 ×{tan−1 [|𝑦|(𝑐2𝑡2−𝑑2)1∕2 𝑧𝑐𝑡 ] + tan−1 [𝑧(𝑐2𝑡2−𝑑2)1∕2 |𝑦|𝑐𝑡 ]} − 3 𝑧 𝑑 𝑐2𝑡2+𝑦2 𝑑2tan−1 [(𝑐2𝑡2−𝑑2)1∕2 𝑧] − 3 |𝑦| 𝑑 𝑐2𝑡2+𝑧2 𝑑2tan−1 [(𝑐2𝑡2−𝑑2)1∕2 |𝑦|] +(2𝑐2𝑡2 𝑑2+ 1)(𝑐2𝑡2 𝑑2− 1)1∕2},(20) and, finally, 𝑄(𝑥, 𝑦, 𝑧, 𝑡) = 𝑐𝑥𝑦H(𝑥)H(𝑦) 4(𝑐𝑡 −𝑧)2H(𝑐𝑡 −𝑧).(21) References [1] Ruehli AE. Inductance calculations in a complex integrated circuit environment. IBM J Res Dev 1972;16(5):470–81. [2] Ruehli AE, Antonini G, Jiang L. Circuit oriented electromagnetic modeling using the PEEC techniques. John Wiley & Sons, Inc., Hoboken, New Jersey; 2017. [3] de Hoop AT. Handbook of radiation and scattering of waves. London, UK: Academic Press; 1995. [4] Antonini G, Orlandi A, Ruehli A. Analytical integration of quasi-static potential integrals on non-orthogonal coplanar quadrilaterals for the PEEC method. IEEE Trans Electromagn Compat 2002;44(2):399–403. [5] Ruehli AE, Antonini G, Esch J, Ekman J, Mayo A, Orlandi A. Non-orthogonal PEEC formulation for time and frequency domain EM and circuit modeling. IEEE Trans Electromagn Compat 2003;45(2):167–76. [6] Lombardi L, Antonini G, Ruehli AE. Analytical evaluation of partial elements using a retarded Taylor series expansion of the Green’s function. IEEE Trans Microw Theory Tech 2018;66(5):2116–27. [7] Wilton D, Rao S, Glisson A, Schaubert D, Al-Bundak O, Butler C. Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains. IEEE Trans Antennas Propag 1984;32(3):276–81. [8] Järvenpää S, Taskinen M, Ylä-Oijala P. Singularity extraction technique for integral equation methods with higher order basis functions on plane triangles and tetrahedra. Internat J Numer Methods Engrg 2003;58(8):1149–65. [9] Järvenpää S, Taskinen M, Ylä-Oijala P. Singularity subtraction technique for highorder polynomial vector basis functions on planar triangles. IEEE Trans Antennas Propag 2006;54(1):42–9. [10] Štumpf M, Antonini G, Ruehli AE. Cagniard-Dehoop technique-based computation of retarded partial coefficients: The coplanar case. IEEE Access 2020;8:148989–96. [11] Štumpf M, Loreto F, Pettanice G, Antonini G. Cagniard–DeHoop technique-based computation of retarded zero-thickness partial elements. Eng Anal Bound Elem 2022;137:56–64. [12] de Hoop AT. A modification of Cagniard’s method for solving seismic pulse problems. Appl Sci Res 1960;B(8):349–56. Engineering Analysis with Boundary Elements 149 (2023) 86–91 91 M. Stumpf et al. [13] de Hoop AT. Large-offset approximations in the modified Cagniard method for computing synthetic seismograms: a survey. Geophys Prospect 1988;36(5):465– 77. [14] de Hoop AT. Reflection and transmission of a transient, elastic, line-source excited SH wave by a planar, elastic bonding surface in a solid. Int J Solids Struct 2002;39(21):5379–91. [15] Štumpf M. Time-domain electromagnetic reciprocity in antenna modeling. Hoboken, NJ: IEEE Press–Wiley; 2019. [16] Štumpf M. Metasurface electromagnetics: The Cagniard-DeHoop time-domain approach. London, UK: IET; 2022. [17] Phillips JR, White JK. A precorrected-FFT method for electrostatic analysis of complicated 3-D structures. IEEE Trans Comput-Aided Des Integr Circuits Syst 1997;16(10):1059–72. [18] Polimeridis AG, Villena JF, Daniel L, White JK. Stable FFT-JVIE solvers for fast analysis of highly inhomogeneous dielectric objects. J Comput Phys 2014;269:280–96. [19] Yucel AC, Georgakis IP, Polimeridis AG, Bağci H, White JK. VoxHenry: FFTaccelerated inductance extraction for voxelized geometries. IEEE Trans Microw Theory Tech 2018;66(4):1723–35. [20] Torchio R, Lucchini F, Schanen J-L, Chadebec O, Meunier G. FFT-PEEC: A fast tool from CAD to power electronics simulations. IEEE Trans Power Electron 2021;37(1):700–13. [21] Lombardi L, Tao Y, Nouri B, Ferranti F, Antonini G, Nakhla MS. Parameterized model order reduction of delayed PEEC circuits. IEEE Trans Electromagn Compat 2019;62(3):859–69. [22] Lombardi L, Loreto F, Ferranti F, Ruehli A, Nakhla MS, Tao Y, et al. Timedomain analysis of retarded partial element equivalent circuit models using numerical inversion of Laplace transform. IEEE Trans Electromagn Compat 2021;63(3):870–9. [23] Loreto F, Romano D, Stumpf M, Ruehli AE, Antonini G. Time-domain computation of full-wave partial inductances based on the modified numerical inversion of Laplace transform method. IEEE Trans Signal Power Integr 2022;1:32–42. [24] Loreto F, Pettanice G, Antonini G, Gad E, Nakhla MS, Tao Y, et al. Modified numerical inversion of Laplace transform methods for the time-domain analysis of retarded partial elements equivalent circuit models. IEEE Trans Electromagn Compat 2021;64(6):2179–88. [25] Lager IE, van Berkel SL. Finite temporal support pulses for EM excitation. IEEE Antennas Wirel Propag Lett 2017;16:1659–62. [26] Štumpf M. Electromagnetic reciprocity in antenna theory. Hoboken, NJ: IEEE Press–Wiley; 2018.