Low Computational Cost Mirror Kirchhoff Approximation for Predicting Shadowing Effect

This paper proposes a 2-dimensional low computational cost mirror Kirchhoff approximation (MKA) and the design of simulation parameters to accurately predict the shadowing gain for an arbitrarily shaped conductor cylinder. The disadvantages of the conventional MKA, such as lacking designs for the simulation parameters and a limited applicable range, have motivated the establishment of an extended MKA. The authors propose the design of the simulation parameters for MKA. The applicable range of MKA is extended to an arbitrarily shaped cylinder by multiple planes. This work finds that only the space domain of the zeroth plane and the angular spectrum domain of the last plane need separate windowing functions for accuracy and the calculation time. The details of those windowing functions are introduced, and their necessities are explained. The authors validate the proposed method for an elliptical conductor cylinder with the size of the human body at millimeter waves (17 GHz – 66.5 GHz). Simulations are conducted by changing the object’s location, direction, and frequencies. The results show that the proposed method presents a good accuracy, which has a low root-mean-square error of less than 0.5 dB by comparing it with a full-wave approach based on the method of moment. Furthermore, the calculation time is improved by 1.4 – 67.2 times by comparing it with the uniform theory of diffraction using special functions.


I. INTRODUCTION
P REDICTION techniques for forward scattering problems are widely used in wireless communication [1] and radar studies [2]. With the widespread use of the 5 th generation mobile communication system (5G), millimeterwave (mmWave) band radio has been widely used [3]. In mmWave, a shadowing effect significantly influences cellphone link performance [4], and hence the prediction technique of the shadowing effect is needed. On the other hand, the enhancement of the forward scatter radar is used to detect low-signature targets in many present works [5]- [7]. It is important to predict forward scattering accurately and quickly for evaluating radar cross sections (RCSs) [6].
As the current prediction methods of forward scattering, there are mainly two types of prediction techniques, i.e., the empirical model and the deterministic model. The empirical models (e.g., log normal shadowing [8], Cheung method [9], and Walfisch-Bertoni model [10]) are fast and easy to use. However, the results of these methods may need to be calibrated by the experimental data, which is not economical. For an economical prediction, the deterministic models are preferable. Conventionally, full-wave electromagnetic approaches (e.g., finite-difference time-domain (FDTD) method [11], method of moment (MoM) [12], multilevel fast multipole algorithm (MLFMA) [13], and finite element method (FEM) [14]), have been used to predict the shadowing effect. These methods solve Maxwell's equations directly and therefore have a good accuracy. However, a high computational cost is a significant drawback for largescale problems. Therefore, another branch of deterministic models, namely, high-frequency asymptotic approximations, is expected to reduce the computational complexity. Highfrequency asymptotic approximation can be divided into raybased and source-based methods.
Ray-based methods, including geometrical optics approximation (GO) [15], geometrical theory of diffraction (GTD) [16], and uniform theory of diffraction (UTD) [17], consider electromagnetic fields as rays. The methods based on Fermat's principle have closed-form analytic solutions, resulting in a lower computational cost. GO describes the propagation phenomena including incidence, reflection, and refraction, but not diffraction. GTD introduces diffraction around edges and smooth objects, while its coefficients are divergent at the shadowing boundary. Its extension UTD uses the Fresnel integral to overcome the discontinuity of GTD near the shadowing boundary. However, when UTD calculates a creeping wave diffraction [18] from an elliptical cylinder, not only the Fresnel integral but also the ellipse integral needs additional calculation time. Furthermore, for a complex shaped object, since the reflection point of the curved surface cannot be easily found by the imaging method, reflection point searching requires plenty of time in the ray-based method.
Source-based methods, such as physical optics approximation (PO) [19] and Kirchhoff approximation (KA) [20], consider electromagnetic fields as equivalent sources. The methods based on Huygens' principle numerically calculate the scattered field for a complex shaped region and hence have a good balance between accuracy and computational complexity. In the PO calculation, the scattered field is calculated from the shadowing object and can be split into reflection radiation to the specular direction and shadow radiation to the forward direction [21]. For the forward scattering problem, the contribution is mainly determined by the shadow radiation. However, since the shadow radiation is calculated by the integration over the cross-section of the geometrical shadow region [22], PO has difficulty evaluating the influence of the thickness on forward scattering problem. In the KA calculation, the scattered field is calculated from the vacuum region around the shadowing object, which also does not consider the thickness of the shadowing object. Although physical theory of diffraction (PTD) [23] combines the advantages of PO and GTD to deal with the forward scattering problem for a thick object, both special functions and numerical integrals require substantial computational costs. For a lower computational cost, the mirror Kirchhoff approximation (MKA) [24] based on KA is used to predict the forward scattering problem for a thick object. Although MKA is the numerical approach, which may be slower than the analytic approach, it is still extremely fast due to the lower computational complexity of the angular spectrum method (ASM) [25]. ASM uses the fast Fourier transform (FFT) [26] to transfer the field between the space and angular spectrum domains.
However, current MKA has not yet designed the parameters of FFT (e.g., interval, windowing function, and windowing size in the space or angular spectrum domains) and is a special case of a thick rectangular object. A small FFT size is expected to reduce the computational cost. Although the formulation processes are quite different, MKA based on FFT is mathematically the same as the wide-angle split-step Model of the angular spectrum method.
parabolic equation (SSPE) [27], [28] or beam propagation method (BPM) [29]. In the SSPE and BPM calculation, Tukey and Hanning windows are widely discussed for forward scattering problems [30]- [34], while they need a large truncation region resulting in a high computational cost.
On the other hand, in the PO calculation, the windowing function based on the Fresnel zone number (FZN) [35]- [37] can achieve a small truncation region within an acceptable accuracy, while those works only consider backward scattering problems. Moreover, none of the above works have designed the windowing function and interval in the angular spectrum domain. These factors affect the computational cost and accuracy. Therefore, a lack of the designed parameters of the FFT is a common issue for all FFT-based methods. This research aims to develop a low computational cost MKA with the designed parameters of FFT for the forward scattering problem. This paper proposes the design of the truncation region and the resolution of FFT. In the space domain, the windowing function based on FZN [35]- [37] is extended to the forward scattering problem. In the angular spectrum domain, a rectangular windowing function is proposed for both accuracy and computational cost. The applicable range of MKA [24] is extended to the arbitrarily shaped cylinder via multiple planes. The paper also proposes a combination of windowing functions in those multiple planes for a better accuracy and lower calculation time. The truncation region, resolution and combination of the windowing functions proposed in this study can also be applicable to SSPE and BPM for fast calculations with an acceptable accuracy.
The remainder of this paper is organized as follows. The design of simulation parameters for MKA is explained in Section II. The application of MKA for an arbitrarily shaped object is shown in Section III. The simulation conditions and methods are described in Section IV. The results are presented and discussed in Section V. Finally, a conclusion of this work is given in Section VI.

II. DESIGN OF SIMULATION PARAMETERS FOR MKA
This section explains the design of the simulation parameters for MKA. A 2D problem is considered by assuming that the problem is uniform along the y-axis. MKA in [24] uses the angular spectrum method (ASM) [25] by applying an inverse Fourier transformation (IFT) as (1) and a Fourier transformation (FT) as (2). Figure 1 shows that ASM calculates the propagation between two parallel planes z = ld and z = (l + 1)d as where x is the parameter of the x-space domain. k x and k z are the parameters of the angular spectrum domain corresponding to the x-space and z-space domains, respectively. E l,i (·) is the total scalar electric field in the l th plane z = ld (l = 0, 1, ..., L) for region i (i = 1, 2). Regions 1 and 2 are the regions above and below the shadowing object, respectively, which are the vacuum regions with the freespace wave number k 0 .Ẽ l,i (·) is the IFT of E l,i . d is the distance between two planes. The cases of |k x | ≤ k 0 and |k x | > k 0 in (3) represent the propagation and evanescent waves, respectively. To apply FFT for a faster calculation speed, discretization for the numerical calculation and truncation for the integral over the finite interval are needed in the space and angular spectrum domains. Suppose the range of FFT and the sampling interval in the space domain are T and ∆x, respectively, and the sampling interval in the angular spectrum domain is 2π/T . Since the FFT size is determined by T /∆x, there are two ways to reduce the calculation complexity, i.e., increasing ∆x and decreasing T .
Normally, [−T /2, T /2] and [−k 0 , k 0 ] are considered as the ranges of truncated regions in the space and angular spectrum domains, respectively. However, by using the Nyquist [38] and λ/10 sampling criteria, the authors propose [−a, a] (2a < T ) and [−k w , k w ] (k w < k 0 ) as the ranges of truncated regions in the space and angular spectrum domains, respectively, for a better accuracy and lower calculation cost.
Section II-A proposes ∆x, which considers the Nyquist sampling criterion for including the evanescent wave. Section II-B proposes a and updates ∆x to ∆x , which considers the λ/10 sampling criterion. Section II-C and Section II-D propose T and k w , respectively.

A. DISCRETIZATION IN THE SPACE DOMAIN
This subsection discusses the spatial resolution. Equation (1) decomposes the source field E l,i as k x -varying plane waves with a complex amplitude ofẼ l,i . Then, those plane waves propagate to the following plane and combine to the aperture field E l+1,i by (2). Evanescent waves exist in those plane waves with |k x | > k 0 since the complex exponential factor e −jkzd becomes the decay exponential factor e − √ k 2 x −k 2 0 d , according to (3). It is known as an evanescent wave where the energy is not transferred but decays exponentially. Usually, the calculation ignores the evanescent wave when d is large [39]. However, when the value of d is small, the evanescent wave with the value of e −   (4), and the Nyquist sampling criterion ∆x ≤ π/ max(k x ) calculates the proposed ∆x by (5).

B. TRUNCATION IN THE SPACE DOMAIN
The proposed ∆x may not satisfy the λ/10 sampling criterion, which is widely used for accuracy [24]. The meaning of the λ/10 sampling criterion is that there are at least ten sampling points in one phase period (one cycle of phase rotation), which is sufficient. This subsection formulates the size of the spatial truncation region and updates the spatial resolution when considering the λ/10 sampling criterion. Figure 2 shows the model of truncation by the spatial window. An object shadows the line-of-sight (LoS) path, which is defined as the z-axis, between a transmitter (Tx) and a receiver (Rx). The point of the surface of the object closest to Tx is defined as the zeroth plane z = 0, which is perpendicular to the z-axis. The distances from the zeroth plane to Tx and Rx are b and c, respectively. w 0,i (i = 1, 2) is the x-coordinate of two object edges in the zeroth plane. a is the window size. The cylindrical wave incident to the zeroth plane can be calculated using the Hankel function of the second kind H the Fresnel region approximation [40] is applied as (6).
Equation (6) gives two key intuitions (P1-P2): The complex amplitude term H  approximately x 0. Truncation by a rectangular windowing function will cause a boundary discontinuity, which takes the extra truncation error in the edge. Thus, there is a need for a decay windowing function to reduce the truncation error.
• P2: Decreasing Phase Period in the Space Domain The phase term e −j k 0 2b x 2 is x 2 -dependent. Thus, the phase period around x 0 is much larger than λ and decreases when x is larger. However, the phase period converges to λ.  [24] to visualize the behavior of the above key intuitions P1-P2. The black line is the plot with the interval of λ/10 as a reference. The red points are the plot with the interval of the proposed ∆x in (5), which is larger than λ/10. We can observe that there are more than ten red sampling points in the phase period near the z-axis. Thus, for the shadowing problem where a region of interest is near the edge, the interval of ∆x near the edge can satisfy ten sampling points within one phase period there.
The region far away from the edge, where the integral of the rapidly oscillating function E 0,i (x)e jkxx does not contribute to the final integral due to the cancellation, can be truncated. Those regions can be numerically integrated in vain if the number of samples per phase period is sufficiently large, e.g., 10 sampling points per phase period. However, the interval of ∆x in the phase period far away from the edge may not provide sufficient sampling points according to key intuition P2. The numerical integration of the vanishing region is erroneous when the number of samples per phase period is insufficiently small. Therefore, it is better to truncate rather than inaccurately integrate the region numerically.
Considering key intuition P1 for avoiding the discontinuity by truncation, the authors choose the windowing function [35] based on the Fresnel zone number by (7).
where W Sp l,i (·) is the spatial windowing function with the edge of w l,i . n(·) is the Fresnel zone number.
The radiation integral can be cancelled if the phase change number n p ≥ 7 (at least a 7π phase change) [35]. Figure 3(b) is the windowed truncation of Fig. 3(a) with the values of w 0,1 = −w 0,2 = 0.25 m and n p = 7. The authors update ∆x of (5) to ∆x as (9), which can satisfy each phase period within the window that should be sampled by at least n s sampling points (n s = 10 for λ/10 sampling criterion). The derivation of (9) is shown in Appendix A.
For the LoS case, w 0,i is replaced by 0 in (7)- (9). Appendix A also derives the spatial windowing size a as (10).

C. DISCRETIZATION IN THE ANGULAR SPECTRUM DOMAIN
This subsection formulates the angular spectrum resolution. Assuming w 0,i , a, x b, c, a Taylor expansion can be used to approximate (11). The windowed field distributed on the region |x − w 0,i | ≤ a of the zeroth plane is calculated by KA as (12).
substituting (12) into (1), the Fresnel integrals S(·), C(·) are used to analytically calculate an asymptotic form ofẼ 0,i (k x ). Due to space, only the real part ofẼ 0,i (k x )/H Although (13) can only provide a rough approximation compared with numerical integration, it still gives three key intuitions (P3-P5): • P3: Amplitude Peaks in the Angular Spectrum Domain From (13), we can determine that k x = k 0 a/b and k x = k 0 w 0,i /b are the poles of the asymptotic form ofẼ 0,i (k x ). Those values correspond to the directions of incident rays to the edges of window and the object in ray tracing. The peaks of oscillatingẼ 0,i (k x ) around those poles make a significant contribution to the integral.
Thus, the phase periods with sizes of 2π/|w 0,i | and 2π/a are constant everywhere of k x .
• P5: Decay Amplitude in the Angular Spectrum Domain The envelope ofẼ 0,i (k x ) is k −1 x -dependent.Ẽ 0,i (k x ) decays fast when k x is far away from peaks. This characteristic influences the selection of the windowing function. The above gives the characteristics ofẼ 0,i (k x ). However, for conducting the next FT by (2), the complex factor E 0,i (k x )e −jkzd should be analyzed. The authors approximate the complex factor e −jkzd around k x 0 as (14).
Equation (14) is identical for the Fresnel region approximation (6), but in the angular spectrum domain. This approximation gives key intuition P6: The phase of e −jkzd is k x 2 -dependent near the line of sight. Thus, e −jkzd has a deceasing phase period with an increase in |k x |. The phase ofẼ 0,i (k x )e −jkzd can be classified into two cases, C1 and C2, by using (13) and (14).
• C2: Decreasing Phase Period in the Angular Spectrum Domain for a Large d In case of a large d, factor e −jkzd dominates the phase ofẼ 0,i (k x )e −jkzd . Thus, by key intuition P6, the phase period is decreasing in the angular spectrum domain.
Using the above phase information, the authors propose T corresponding to the angular spectrum resolution as follows.
From key intuition P3, the phase period at the poles should be sampled by n s points for accuracy in the angular spectrum domain. For the LoS condition, k x = 0 corresponding to the directions of the direct ray is the main peak. For the NLoS condition, the poles ofẼ 0,i (k x )e −jkzd are k x = k 0 a/b and k x = k 0 w 0,i /b, which are the same as the poles of E 0,i (k x ) due to |Ẽ 0,i (k x )e −jkzd | = |Ẽ 0,i (k x )|. Additionally, for case without spatial windowing truncation, since the terms including a in (13) vanish for a → ∞, the pole of k x = k 0 a/b and the phase term of ak x disappear. Using the phase information of cases C1 and C2, the authors propose T as (16).
where T 1 is the proposed T in case C1 or case C2 for the LoS condition with spatial windowing truncation. T 2 is the proposed T in case C2 for the NLoS condition with spatial windowing truncation. T 3 is the proposed T in case C1 or case C2 for the LoS condition without the spatial windowing truncation. T 4 is the proposed T in case C2 for the NLoS condition without the spatial windowing truncation. The derivation of (16) is shown in Appendix B.

D. TRUNCATION IN THE ANGULAR SPECTRUM DOMAIN
This subsection formulates the size of the angular spectrum truncation region. In case C1, there is no need to consider the extra discretization error by the equally sampled points, since the phase period is almost constant. However, there is a vast discretization error for the equally sampled points in case C2 since the phase period decreases with an increasing |k x |. The region far away from the poles in the angular spectrum domain, where the integral of the oscillating function does not contribute to the final integral, should be truncated rather than inaccurately integrate the insufficiently sampled region numerically. The angular spectrum domain is truncated by using the rectangular windowing function for case C2. The rectangular windowing function is chosen because the amplitude ofẼ 0,i (k x )e −jkzd is |Ẽ 0,i (k x )|, which decays fast according to key intuition P5. Thus, a rectangular windowing function results in a negligibly small truncation error. Figures 4(a) and (b) show the examples before and after the truncation with the value of d = 7.85 m in the same condition as Fig. 3(b) [24]. The windowing function for the angular spectrum domain is (17).
where W An (·) is the angular spectrum windowing function. k w is the windowing size in the angular spectrum domain and is proposed as (18).
where n c is the Nyquist sampling number per phase period, i.e., 2. The derivation of (18) is shown in Appendix B. Appendix B also explains the reason for using the Nyquist sampling criterion in the angular spectrum domain, which is different from the λ/10 sampling criterion in the space domain.

III. APPLICATION OF MKA FOR AN ARBITRARILY SHAPED OBJECT.
This section introduces the application of MKA for an arbitrarily shaped object. The arbitrarily shaped cylinder is approximated to the combination of several rectangular cylinders, as shown in Fig. 5. These rectangular cylinders can be seen as slices of an arbitrarily shaped cylinder by multiple planes, as shown in Fig. 6. In Fig. 6, the distances b and c−Ld are large, while d is small. There are (L + 1) planes, and all the planes need to be perpendicular to the Tx -Rx line [41]. The interval between two planes is determined by the maximum propagation angle θ m [27] as (19).
By repeatedly applying MKA among those planes as (20), the scattered fields can be calculated to evaluate the shadowing gain.
where E MKA l+1,i (·) is the field calculated by MKA [24] in the (l + 1) th plane. w l+1,i (i = 1, 2) are the coordinates of two edges in the (l + 1) th plane. The second term on the righthand side is the mirror image. The fields E MKA l+1,i (·) for the regions i = 1, 2 are calculated separately, until they reach the last plane where two fields are merged.
The new finding is that, only the space domain of the zeroth plane and the angular spectrum domain of the last plane need their respective windowing functions, while other planes or domains do not. To prove this finding, there is a need to clarify what condition needs or does not need the windowing function. According to Section II, the meaning of the windowing function is to truncate the region where they cannot contribute to the integral. However, those regions have discretization errors caused by the equal sampling interval in a decreasing phase period. Thus, the constant phase period does not need a windowing function, and vice versa for the decreasing phase period. The authors explain the characteristics of the phase period for each plane in both domains as follows.
For the zeroth plane, the Hankel function in the space domain has a decreasing phase period according to key intuition P2. Thus, there is a need to take the spatial windowing function. However, since the change of the phase period in the small truncation region is insignificant, the phase period of the windowed Hankel function E 0,i (x) is constant in the truncation region. On the other hand, in the angular spectrum domain,Ẽ 0,i (k x )e −jkzd has a constant phase period according to key intuition P4. Since d is small,Ẽ 0,i (k x ) has a constant phase period according to case C1. Thus, there is no need to take the angular spectrum windowing function.
For the first plane, sinceẼ 0,i (k x ) dominatesẼ 0,i (k x )e −jkzd according to case C1, E 1,i (x) has almost the same constant phase period with E 0,i (x) as (21). Thus, there is no need to take the spatial windowing function. The same approach is used with (13) to analytically derive the above characteristics as (22). Equation (22) shows that E 1,i (x) has a constant phase period.
Considering (21) as the initial condition and repeating (23) and (24) inductively, the authors prove that E l+1,i (x) andẼ l+1,i (k x )e −jkzd have constant phase periods, such as E 0,i (x) andẼ 0,i (k x )e −jkzd , respectively. Thus, there is no need for the windowing function in the space or the angular spectrum domain until the last plane.
Finally, for the last plane in the angular spectrum domain, since e −jkz(c−Ld) dominatesẼ l,i (k x )e −jkz(c−Ld) for the large c − Ld,Ẽ l,i (k x )e −jkz(c−Ld) has a decreasing phase period, according to case C2. Thus, there is a need to take the angular spectrum windowing function.
Above all, we can conclude that only the space domain of the zeroth plane and the angular spectrum domain of the last plane need their separate windowing functions in case C2. Figure 7 shows the flow chart of MKA calculation. The green block corresponds to the proposed formulation men- tioned in Section II. The red remarks are the findings mentioned in Section III. The order of the calculation time for MKA with the proposed parameters is (25).
where N = T α /∆x . t is the calculation time. L is the total number of approximated rectangles. ∆x is calculated by (9). T α is calculated by (16).

IV. SIMULATION
The authors validated the proposed method for an elliptical conductor cylinder, which models the cross section of the human body. As shown in Fig. 8, an elliptical conductor cylinder is placed between Tx and Rx. An electric line source with the cylindrical wave in the perpendicular polarization is considered at mmWave. r 1 and r 2 denote the semi-major and semi-minor axes of the ellipse, respectively. f denotes the frequency. d 1 denotes the horizontal distance between the center of the ellipse and Tx. d 2 denotes the horizontal distance between the center of the ellipse and Rx. ∆d denotes the offset between the center of the ellipse and the Tx-Rx line. θ denotes the rotation angle of the ellipse. There are two simulation scenarios. Scenario 1 varies ∆d from −100λ to 100λ with an interval of 0.5λ. Scenario 2 varies f from 17 GHz to 66.5 GHz with an interval of 0.5 GHz at ∆d = 0 m. Considering the human-size shadowing problem, the authors set the parameters of each scenario as shown in Table 1 and Table 2, respectively. For each scenario, three rotation angles are considered, i.e., θ = 0 • , θ = 45 • , and θ = 90 • .
The processor of the calculating computer is an Intel(R) Core(TM) i7-8750H CPU @ 2.20 GHz. The usable installed memory of the calculating computer is 15.8 GB. The system type of the calculating computer is a 64-bit operating system with a x64-based processor. The simulation software is MATLAB.
Computation time of MKA is compared with that of UTD. To discuss the calculation time later, the authors give the detailed UTD simulation approaches in Appendix C.
MoM using the electric field integral equation as a reference of accuracy is simulated. The piecewise constant function is selected as the basis function of the MoM, and the point matching method is implemented. The mesh size of MoM is set to λ/10 for accuracy.

V. RESULTS AND DISCUSSIONS
This section shows the simulation results, including the shadowing gain and the calculation time in Section V-A. A discussion for evaluating the results is provided in Section V-B.

A. RESULTS
For scenario 1, the shadowing gain results of UTD, MKA and MoM are shown in Fig. 9. The horizontal axis is ∆d (λ). The vertical axis is the shadowing gain results on a decibel (dB) scale. The results show that the accuracy of the proposal is validated by a comparison with MoM for varying the location of the object.
The relative calculation time of MKA is shown in Fig. 10. The horizontal axis is ∆d (λ). The vertical axis is the relative calculation time referenced by UTD. The improvement in the calculation time is compared with UTD for varying the location of the object.
For scenario 2, the shadowing gain results of UTD, MKA and MoM are shown in Fig. 11. The horizontal axis is Plots of the shadowing gain for scenario 1. f (GHz). The vertical axis is the shadowing gain results on a dB scale. The results show that the accuracy of the proposal is validated by a comparison with MoM for varying the frequencies.
The relative calculation time of MKA is shown in Fig. 12. The horizontal axis is f (GHz). The vertical axis is the relative calculation time. The improvement in the calculation time is compared with UTD for varying the frequencies.
Considering MoM as the reference, the authors calculate the root-mean-square error (RMSE) by (26).
where SG MoM j is the shadowing gain calculated by MoM on a dB scale for the j th test and SG Method j is the shadowing gain calculated by UTD or MKA on a dB scale for the j th test. m is the total number of tests.
The comparisons of the RMSE and the absolute computational time between UTD and MKA for each scenario are shown in Table 3 and Table 4, respectively.

B. DISCUSSIONS
The results showed that the proposed low computational cost MKA achieves a good accuracy with a low root-mean-square error of less than 0.5 dB, compared with the MoM as a reference.
The authors evaluated the dominant parts of the calculation time for UTD and MKA as follows. For the LoS case of f = 66.5 GHz, θ = 0 • , ∆d = 50λ, there was one reflection 10   and reflection point searching. Furthermore, the authors evaluated the absolute calculation time for the double diffraction mentioned in [24]. Comparing the calculation time with 210.0 ms for one test in [24], the low computational cost MKA provided a faster calculation time with 5.0 ms for one test.

VI. CONCLUSION
This paper proposed a design for the MKA simulation parameters. The applicable range of MKA was extended to an arbitrarily shaped cylinder through multiple planes. This work found that, only the space domain of the zeroth plane and the angular spectrum domain of the last plane needed their separate windowing functions for both accuracy and the calculation time. The authors validated the proposed method for an elliptical conductor cylinder with the size of the human body at the mmWave (17 GHz -66.5 GHz). The results showed that the proposed method presented a good accuracy with a low RMSE of less than 0.5 dB, compared with the MoM as the reference. Furthermore, the calculation time was improved by 1.4 -67.2 times compared with the UTD using special functions. The designed simulation parameters and the combination of windowing functions presented in this paper can also be applicable to SSPE or BPM for a better accuracy and lower computational cost. Extensions of the applicable range to other opaque object, such as the human body, and 3-dimensional object are future topics. .

APPENDIX A DERIVATIONS OF THE RESOLUTION AND WINDOWING SIZE IN THE SPACE DOMAIN
According to key intuition P2, the phase function φ 1 (·) of the Hankel function around x 0 is The spatial windowing boundary with the coordinates of a i should be n p phase rotations away from the edge with the coordinates of w 0,i .
The coordinate a i , which is the 2π phase smaller than a i , is Thus, the size ∆a i of the last phase period within the window is Since ∆a i should be sampled by n s points, the sampling interval ∆x i should be The proposed interval ∆x should satisfy both ∆x of (5) and ∆x i as Since the updated ∆x may be smaller than ∆x i , for a better accuracy, the window size with sufficient sampling points can be larger than a i , which just satisfies the minimum requirements of truncation with the phase change number of n p . Suppose a is the updated coordinate of the windowing boundary. The coordinate a , which is the 2π phase smaller than a, is Thus, the size ∆a of the last phase period within the window is Since ∆a should be sampled by n s points, a can be derived as

APPENDIX B DERIVATIONS OF THE RESOLUTION AND WINDOWING SIZE IN THE ANGULAR SPECTRUM DOMAIN
In case C1, according to key intuition P4, the size of the phase period everywhere for k x is 2π/|w 0,i | and 2π/a (a > w 0,i ). Each phase period should be sampled by at least n s points as If there is no spatial windowing function, since the term a vanishes, the size of the phase period is 2π/|w 0,i |. In case C2, according to key intuition P6, the phase function φ 2 (·) of e −jkzd around k x 0 is For the NLoS condition in case C2, the phase period at poles k x = k 0 a/b and k x = k 0 |w 0,i |/b should be sampled by n s points.
If there is no spatial windowing function, then the pole of k x = k 0 a/b vanishes.
For the LoS condition in case C2, the phase period of E 0,i (k x )e −jkzd at peak k x = 0 is the same asẼ 0,i (k x ) due toẼ 0,i (k x )e −jkzd | kx=0 =Ẽ 0,i (k x )| kx=0 . Thus, the LoS condition in case C2 has the same conclusion as case C1. Truncation size k w is calculated by using the Nyquist sampling number n c as The authors explain the reason for using the Nyquist sampling criterion in the angular spectrum domain. Figure 13 uses poles k x = ±k 0 a/b and boundaries k x = ±k w to separate Fig. 4(b) into five regions, i.e., A, B, C, D, and E. Regions A and E cannot be sampled by at least n c points. Regions B and D are sampled from n c to n s points. Region C is sampled by at least n s points. The insufficiently sampled regions A and E can be truncated due to the cancellation of the radiation integration. The sufficiently sampled regions C should remain due to the contribution of the integral. However, transition regions B and D, where many peaks of the oscillating 2 i=1Ẽ 0,i (k x )e −jkzd are around the poles, also play an important role in the integration. If we suddenly truncate by using the n s sampling criterion, then transition regions B and D will vanish causing the accuracy issue due to the lack of peaks in |k x | ≥ k 0 a/b. The reason for using the n c sampling criterion for truncation is because it can retain many peaks with heights over n c /n s of the highest peak. The proof is as follows.
Suppose the peak at pole k x = k 0 a/b is the highest peak with a height of h 1 , and the height of the peak at boundary k x = k w is h 2 . According to key intuition P5, since the envelope ofẼ 0,i (k x )e −jkzd is k −1 x -dependent, we have In addition, the sampling numbers of the phase periods at k x = k 0 a/b and k x = k w are n s and n c , respectively. Thus, we have From (42) and (43), we can prove that For the values of n c = 2 and n s = 10, only the peaks with heights less than 20% of h 1 , which do not significantly influence the next FT, are truncated. Thus, the accuracy of the truncation is not an issue. For discretization, it is better not to unify the sampling criteria in the space and angular spectrum domains. If we unify n c = 10 in the angular spectrum domain and want to retain the same peaks with heights over 20% of h 1 , then n s will reach 50 points. According to (16), T α will be 5 times larger than before, causing a significant increase in the calculation time. If we unify n s = 2 in the space domain, according to (10), where the dominant part is a = n 2 s k 0 ∆x 2 + 4πb 2n s k 0 ∆x ∼ 2πb n s k 0 ∆x (45) a will be 5 times larger than before. According to (16b), T 2 will also be 5 times larger than before, causing a significant increase in the calculation time. This paper has already shown that the combination of n c = 2 and n s = 10 can present a good accuracy with a low RMSE of less than 0.5 dB. Therefore, there is no need to unify the sampling criteria in the two domains for a lower computational cost.