Electromagnetic simulations of large and complex targets are frequently addressed using high-frequency algorithms, notably in scenarios such as ship-sea complexes, aircraft, and inlet targets. Among the efficient algorithms for solving cavity scattering problems, the iterative physical optics (IPO) method iteratively refines the induced current of the physical optics (PO) method, enhancing the accuracy in targets with multiple scattering mechanisms. Furthermore, compared with the full-wave algorithm’s discrete size of (100 to 144) facets$ /\lambda^2 $, the IPO method only requires a density of (9 to 25) facets$ /\lambda^2 $ for discretization. This enables it to achieve a larger simulated size-to-wavelength ratio, where λ represents the wavelength. The multilevel fast multipole algorithm (MLFMA) [1] effectively reduces the computational complexity of the moment method from $ O\left( {{N^3}} \right) $ to $ O\left( {{N^2}\; {\mathrm{log}} \;{{N}}} \right) $, where $ N $ represents the number of patches. It is employed in solving the electromagnetic problem of three-dimensional complex targets under high-performance server conditions. The unknown quantity associated with this method can surpass billions in simulations involving thousands of wavelength scenarios. Consequently, it is frequently paired with central processing unit (CPU) or graphics processing unit (GPU) parallel technology to address the challenge of low simulation efficiency [2,3]. The PO method exclusively considers the induced current in the lit region and is suitable for the scattering problem of smooth electrically large targets [4,5]. To mitigate the computational complexity of the shadow judgment process, techniques such as graphics processing [6,7] or rendering methods are employed.

On the contrary, utilizing the shotting and bouncing rays (SBR) method [8] based on the discrete size of the aperture surface being a point encounters a significant challenge—the immense number of ray tubes poses a memory limit issue in large scenes. Nevertheless, given the independent tracking processes of each ray tube, addressing the phase of multiple reflections between various patches can resolve the challenge of scenes involving multiple scattering [9-12]. As a solution, GPU parallel tracking is employed to tackle the scattering problem in large aircraft scenarios [11,13-15]. However, after the ray bounces multiple times inside the cavity, the accuracy decreases due to the divergence factor and the divergence conditions of the ray tube. However, the accuracy diminishes after multiple bounces inside the cavity due to the divergence factor and divergence area of the ray tube [8,16]. The IPO method grounded in the iterative calculation of induced current and based on the magnetic field integral equation (MFIE) demonstrates the high simulation accuracy in cavity scattering problems. It has been successfully addressed in both two-dimensional [17] and three-dimensional [18] scattering problems. Moreover, the cost of the number of patches in the IPO method is obviously reduced compared with the full-wave algorithm [2]. The number of iterations required by the IPO method increases significantly as the size of the cavity expands [19]. The forward-backward IPO algorithm effectively solves the convergence challenge associated with multiple iterations of cavity induced current [20].

Considering the visibility function between patches can enhance the iterative convergence of the IPO method. However, for complex targets, relying solely on judging the shadow relationship based on the normal vector ignores the occlusion problem of any shape between the test point and the integration point, inevitably leading to errors. As a consequence, this method proved to be insufficiently accurate for more intricate target scattering calculations [19,21]. Achieving the higher scattering accuracy can be obtained by considering the iterative or superimposed diffraction mechanism contribution [22] between physical theory of diffraction (PTD) current and surface current [23]. The shadow area radiation scheme can reduce the costly computing and storage requirements of shadow judgment, thereby enhancing the convergence speed of the IPO method [24]. For the scattering problem involving arbitrarily shaped cavities, employing a shadow judgment method under strict judgment criteria obtains an accurate visibility function, which can reduce misjudgment in patch current iterations and improve the convergence of the algorithm. However, the complexity of brute traversal judgment between patches will increase the computational complexity from $ O\left( {{N^2}} \right) $ to $ O\left( {{N^3}} \right) $. Shadow judgment based on accelerated tree structures such as octree [25,26] and K-dimensional tree (KD-Tree) [27,28] is applied for verification of ray tracing in the PO and SBR methods. This structure proves beneficial in reducing the need for a high number of iterative judgments, hereby improving simulation efficiency.

In this context, the time harmonic factor of the electromagnetic (EM) field is chosen to be $ {{e}}^{-\mathrm{i}\omega t} $. We consider the establishment process of a bounding volume hierarchy (BVH) based on flat triangle patches. Secondly, the visibility function is taken into account during the cavity scattering process of the tree-shaped acceleration structure. Finally, criteria for judging shadow occlusion between patches are provided. Numerical simulations demonstrate that the proposed method ensures the accuracy of the shadow judgment function between patches and reduces computational complexity from $ O\left( {{N^3}} \right) $ to $ O\left( {{N^2}\;{\mathrm{log}}\; {{N}}} \right) $ through the BVH technology, which is beneficial for expanding the simulation size and wavelength ratio while ensuring the accuracy.

Utilizing the Morton code to construct a binary tree with hierarchical bounding boxes in Fig. 1, the process of shadow occlusion between patches involves initial judgment with the bounding box, followed by examination within the internal patches. Classifying the assessment of the lit-shadow relationship between patches can be framed as a ray occlusion problem. The origin of the ray is defined to $ {\mathbf{P}}\left( {{p_x}{\rm{,}}{\text{ }}{p_y}{\rm{,}}{\text{ }}{p_z}} \right) $, and the nearest intersection point with the box is denoted as $ \mathbf{Q}\left(q_x\rm{,}\text{ }q_y\rm{,}\text{ }q_{\textit{z}}\right) $. Any point $ {{\bf{W}}}({w_x}{\rm{,}}{\text{ }}{w_y}{\rm{,}}{\text{ }}{w_z}) $ along the ray can be expressed as $ {\bf W} = {\mathbf{P}} + t\widehat {{\mathbf{P}}{{\bf{Q}}}} $, where $ t $ is the length of the ray, and $ \widehat {{\mathbf{OP}}} $ represents the unit propagation direction of the ray. Based on the upper bound $ {\bf{R}_{{\mathrm{max}} }}\left( {{\beta _1}{\rm{,}}{\text{ }}{\beta _2}{\rm{,}}{\text{ }}{\beta _3}} \right) $ and the lower bound $ {\bf{R}_{{\mathrm{min}} }}\left( {{\alpha _1}{\rm{,}}{\text{ }}{\alpha _2}{\rm{,}}{\text{ }}{\alpha _3}} \right) $ of the box, we can get

$ {{{\bf{R}}}_{{\mathrm{min}} }}\left( {{\alpha _1}{\rm{,}}{\text{ }}{\alpha _2}{\rm{,}}{\text{ }}{\alpha _3}} \right) \leq {{\bf{W}}} \leq {{{\bf{R}}}_{{\mathrm{max}} }}\left( {{\beta _1}{\rm{,}}{\text{ }}{\beta _2}{\rm{,}}{\text{ }}{\beta _3}} \right).{\text{ }} $ (1)

When the ray intersects the bounding box, $ {t_{{\mathrm{max}} }} $ and $ {t_{{\mathrm{min}} }} $ satisfy

$ {t_{{\rm{min}} }} = {\rm{max}} \left\{ {\frac{{{x_{{\rm{min}} }} - {{\mathbf{Q}}_x}}}{{{{\widehat {{\mathbf{QP}}}}_x}}}{\rm{,}}{\text{ }}\frac{{{y_{{\rm{min}} }} - {{\mathbf{Q}}_y}}}{{{{\widehat {{\mathbf{QP}}}}_y}}}{\rm{,}}{\text{ }}\frac{{{z_{{\rm{min}} }} - {{\mathbf{Q}}_z}}}{{{{\widehat {{\mathbf{QP}}}}_z}}}} \right\} \tag{2a}$ ()

$ {t_{{\rm{max}} }} = {\rm{min}} \left\{ {\frac{{{x_{{\rm{max}} }} - {{\mathbf{Q}}_x}}}{{{{\widehat {{\mathbf{QP}}}}_x}}}{\rm{,}}{\text{ }}\frac{{{y_{{\rm{max}} }} - {{\mathbf{Q}}_y}}}{{{{\widehat {{\mathbf{QP}}}}_y}}}{\rm{,}}{\text{ }}\frac{{{z_{{\rm{max}} }} - {{\mathbf{Q}}_z}}}{{{{\widehat {{\mathbf{QP}}}}_z}}}} \right\}{\mathrm{.}} \tag{2b}$ ()

When the ray intersects the bounding box, it must satisfy $ {t_{{{\mathrm{max}}}}} \ge {t_{{{\mathrm{min}}}}} $ and $ {t_{{{\mathrm{max}}}}} \ge 0 $. Additionally, tracing between rays, the recursive judgment is based on the established binary BVH structure in Fig. 2. First, the ray undergoes an occlusion judgment with the outermost bounding box. If the judgment result indicates an intersection, the process recursively continues with the left subtree bounding box in the next layer until reaching the bottom layer patch that intersects. This process leverages the advantages of the tree structure and avoids the process of judging intersecting patches through traversal in the traditional method in Fig. 3. Consequently, the complexity of the process is reduced from $ O({N^3}) $ to $ O({N^2}\;{\mathrm{log}}\; N) $.

For the cavity scattering problem of a plane wave with an incident direction $ {\widehat {\mathbf{k}}_i} $ in free space, the scattering electric field $ {\mathbf{E}}_{{{\mathrm{all}} }}^s = {\mathbf{E}}_{{{\mathrm{ext}}}}^s + {\mathbf{E}}_{{{\mathrm{int}}}}^s $ originates from the external part of the cavity $ {\mathbf{E}}_{{{\mathrm{ext}}}}^s $ and the internal part of the cavity $ {\mathbf{E}}_{{{\mathrm{int}}}}^s $. The IPO method is employed as the solver for calculating the internal scattered field $ {\mathbf{E}}_{{{\mathrm{int}}}}^s $. A virtual aperture $ {S_a} $ is placed at the open-ended cavity to ensure equivalence of the uniform incident electric field $ {{\mathbf{E}}_i}({\mathbf{r}}) = {{\hat {\bf{E}}}_i}{e^{ - {\mathrm{i}}k{{\widehat {\mathbf{k}}}_i} \cdot {\mathbf{r}}}} $ and the initial equivalent electric and magnetic surface currents $ {\mathbf{J}}_a^0 $ and $ {\mathbf{M}}_a^0 $ at the aperture are obtained from the Kirchhoff approximation as

$ {\mathbf{J}}_a^0({\mathbf{r}}) = {\hat n}\left( {{{\mathbf{r}}_a}} \right) \times {{\mathbf{H}}_i}\left( {{{\mathbf{r}}_a}} \right) = {H_0}{\hat {\bf{n}}} \times {{{\hat {\bf{H}}}}_i} \cdot {{\boldsymbol{e}}^{{\mathrm{i}}k{{\widehat {\mathbf{k}}}_i} \cdot {{\mathbf{r}}_a}}} \tag{3a}$ ()

$ {\mathbf{M}}_a^{0}({\mathbf{r}}) = {{\mathbf{E}}_i}\left( {{{\mathbf{r}}_a}} \right) \times {\hat {\bf{n}}}\left( {{{\mathbf{r}}_a}} \right) = {E_0}{{{\hat {\bf{H}}}}_i} \times {\hat {\bf{n}}} \cdot {{\boldsymbol{e}}^{{\mathrm{i}}k{{\widehat {\mathbf{k}}}_i} \cdot {{\mathbf{r}}_a}}} \tag{3b} $ ()

where $ {\hat {\bf{n}}} $ represents the unit normal vector on the aperture $ {S_a} $ pointing to the inside of the cavity. Considering the occlusion relationship between the aperture patch and the inner wall patch in Fig. 3. A visibility function $ V\left( {{\mathbf{r}}{\rm{,}}{\text{ }}{\mathbf{r}'}} \right) $ is introduced. Specifically, $ V_{{\mathrm{ac}}}\left(\mathbf{r}\rm{,}\text{ }\mathbf{r}'\right)=1 $ represents only if there is no occlusion between the field point $ {\mathbf{r}} $ on the aperture ${S_a}$ and the source point $ {\mathbf{r}'} $ on the cavity wall ${S_w}$, otherwise it is $ V_{\mathrm{ac}}\left(\mathbf{r}\rm{,}\text{ }\mathbf{r}'\right)=0 $. The subscript $a$ represents the patch on the aperture, while the subscript $c$ represents the patch on the inner wall of the open-ended cavity. The initial magnetic field $ {\mathbf{H}}_w^{0}({\mathbf{r}}) $ on the inner wall $ {S_w} $ of the cavity is approximated by the Kirchhoff approximation principle to $ {\mathbf{J}}_a^{0} $ and $ {\mathbf{M}}_a^{0} $ passing through the aperture $ {S_a} $ as follows [18]:

$ $$\begin{array}{c}{\mathbf{H}}_w^0\left({\mathbf{r}}_c\right)={\int }_{S_a}{V}_{\mathrm{ac}}\left({\mathbf{r}}_a,{\mathbf{r}}_c\right){\mathbf{J}}_a^0\left({\mathbf{r}}_a\right)\times \mathrm{\nabla }{G}_0\left({\mathbf{r}}_a,{\mathbf{r}}_c\right)\mathrm{d}{S}_a+\frac{1}{\mathrm{i}k{Z}_0}\mathrm{\nabla }\times {\int }_{S_a}{V}_{\mathrm{ac}}\left({\mathbf{r}}_a,{\mathbf{r}}_c\right){\mathbf{M}}_a^0\left({\mathbf{r}}_a\right)\times \mathrm{\nabla }{G}_0\left({\mathbf{r}}_a,{\mathbf{r}}_c\right)\mathrm{d}{S}_a.\end{array}$$ $ (4)

Here, a strict shadow relationship is considered between the aperture and the inner wall with $ \mathbf{H}_{S_a}^{\mathrm{sha}}\left(\mathbf{r}_c\right)\equiv0 $. The wave impedance of a plane wave in free space is represented as ${Z_0}$. Moreover, $ \nabla {G_0}\left( {{\mathbf{r}}{\rm{,}}{\text{ }}{\mathbf{r}'}} \right) $ represents the gradient of the free-space dyadic Green’s function as follows:

$ \nabla {G_0}\left( {{\mathbf{r}}{\rm{,}}{\text{ }}{\mathbf{r}'}} \right) = \left( {\frac{1}{R} - {\mathrm{i}}k} \right)\frac{{{e^{{\mathrm{i}}kR}}}}{{4\pi R}}{\hat {\bf{R}}} $ (5)

where ${{\bf{R}}}$ represents the position vector between the field point $ {\mathbf{r}} $ and the source point $ {\mathbf{r}'} $. Its unit vector is expressed as$ {\hat {\bf{R}}} = \left( {{\mathbf{r}} - {\mathbf{r}'}} \right)/\left| {{\mathbf{r}} - {\mathbf{r}'}} \right| $, and $ R $ represents the magnitude of the position vector $ {{\bf{r}}} - {{\bf{r}}'} $. In addition, the perfect electrical conductor (PEC) is considered, so the electric field on the inner wall $ {{\mathbf{E}}_w} \equiv {\mathbf{0}} $ and the induced magnetic current $ {{\mathbf{M}}_w} \equiv 0 $. In addition, based on the PO method induced current approximation, the iterative initial value of the IPO method is expressed as $ {\mathbf{J}}_w^{0}\left( {{{\mathbf{r}}_c}} \right) = 2{\hat {\bf{n}}} \times {\mathbf{H}}_w^{0}\left( {{{\mathbf{r}}_c}} \right) $. Generally, the Nth iteration expression of the IPO method can be expressed as follows using MFIE:

$ {\mathbf{J}}_w^N\left( {{{\mathbf{r}}_c}} \right) = {\mathbf{J}}_w^{0}\left( {{{\mathbf{r}}_c}} \right) + 2{\hat {\bf{n}}} \times \int_{{S_{c}}} {{V_{{\mathrm{cc}}}}} \left( {{{\mathbf{r}}_c}{\rm{,}}{\text{ }}{{{\mathbf{r}'}}_c}} \right){\mathbf{J}}_w^{N{ - 1}}\left( {{\mathbf{r}'}} \right) \times \nabla {G_0}\left( {{{\mathbf{r}}_c}{\rm{,}}{\text{ }}{\mathbf{r}'}} \right){{\mathrm{d}}}{S_c}{\mathrm{.}}{\text{ }} $ (6)

The visibility function $ {V_{{\mathrm{cc}}}}\left( {{\mathbf{r}}{\rm{,}}{\text{ }}{\mathbf{r}'}} \right) $ between inner wall patches is introduced here to enhance the accuracy of the induced current iteration. When the inner wall induced current $ {{\mathbf{J}}_N}\left( {{{\mathbf{r}}_c}} \right) $ falls below the convergence threshold $ \varepsilon $, it indicates that the induced current has converged. Typically, the threshold for current convergence $\varepsilon $ is set to 0.01. At this point, the Kirchhoff approximation principle is reintroduced to update the electromagnetic current $ {\mathbf{J}}\left( {{{\mathbf{r}}_a}} \right) $ and $ {\mathbf{M}}\left( {{{\mathbf{r}}_a}} \right) $ at the aperture $ {S_a} $. In this scenario, when the magnetic current $ {\mathbf{M}} $ is constant at 0, the equivalent process can be simplified to

$ {\mathbf{H}}\left( {{{\mathbf{r}}_a}} \right) = {\mathbf{H}}_{{S_c}}^{{{\mathrm{lit}}}}\left( {{{\mathbf{r}}_a}} \right) = \int_{{S_{c}}} {{V_{{\mathrm{ca}}}}} \left( {{{\mathbf{r}}_c}{\rm{,}}{\text{ }}{{\mathbf{r}}_a}} \right){{\mathbf{J}}_N}\left( {{{\mathbf{r}}_c}} \right) \times \nabla {G_0}\left( {{{\mathbf{r}}_a}{\rm{,}}{\text{ }}{{\mathbf{r}}_c}} \right){{\mathrm{d}}} {S_{c}}{\mathrm{,}} \tag{7a}$ ()

$ {\mathbf{E}}\left( {{{\mathbf{r}}_a}} \right) = {\mathbf{E}}_{{{S}_c}}^{{{\mathrm{lit}}}}\left( {{{\mathbf{r}}_a}} \right) = \frac{1}{{{\mathrm{i}}k{Y_0}}}\nabla \times \int_{{S_{{\mathrm{c}}}}} {{V_{{\mathrm{ca}}}}} \left( {{{\mathbf{r}}_c}{\rm{,}}{\text{ }}{{\mathbf{r}}_a}} \right){{\mathbf{J}}_N}\left( {{{\mathbf{r}}_c}} \right) \times \nabla {G_0}\left( {{{\mathbf{r}}_a}{\rm{,}}{\text{ }}{{\mathbf{r}}_c}} \right){{\mathrm{d}}} {S_c}{\mathrm{.}} \tag{7b}$ ()

The visibility function $ {V_{{\mathrm{ca}}}} $ represents the lit-shadow relationship from any inner wall patch to the aperture patch, which is numerically equal to $ {V_{{\mathrm{ac}}}} $ in (4). Finally, the far-field electric field scattering is solved with the Kirchhoff approximation principle with the $ {\mathbf{J}}\left( {{{\mathbf{r}}_a}} \right) = - {\hat {\bf{n}}} \times {\mathbf{H}}\left( {{{\mathbf{r}}_a}} \right) $ and $ {\mathbf{M}}\left( {{{\mathbf{r}}_a}} \right) = - {\mathbf{E}}\left( {{{\mathbf{r}}_a}} \right) \times {\hat {\bf{n}}} $ at the updated aperture when the distance $ r \to \infty $. The admittance $ {Y_0} $ is obtained by $ {Y_0} = 1/{Z_0} $. The far-field scattered electric field $ {{\mathbf{E}}_s} $ can be expressed as

$ $$\begin{array}{c}{\mathbf{E}}_s(\mathbf{r})={\mathbf{E}}_{S_a}(\mathbf{r})=-{\int }_{S_a}\mathbf{M}\left({\mathbf{r}}_a\right)\times \mathrm{\nabla }{G}_0\left(\mathbf{r},{\mathbf{r}}_a\right)\mathrm{d}{S}_a+\frac{1}{\mathrm{i}k{Y}_0}\mathrm{\nabla }\times {\int }_{S_a}\mathbf{J}\left({\mathbf{r}}_a\right)\times \mathrm{\nabla }{G}_0\left(\mathbf{r},{\mathbf{r}}_a\right)\mathrm{d}{S}_a.\end{array}$$ $ (8)

In summary, we have introduced three types of visibility functions $ {V_{{\mathrm{ac}}}} $, $ {V_{{\mathrm{cc}}}} $, and $ {V_{{\mathrm{ca}}}} $, which contribute to enhancing the accuracy of the induced electromagnetic current in the equivalent process. The method avoids the erroneous current iteration steps in the presence of occlusion between patches and facilitate convergence during the current iteration process. The occlusion judgment between patches is detailed in Section 4.

In fact, the results of shadow region judgments are frequently utilized in the equivalent calculations of the aperture surface and the electromagnetic current iteration of the inner wall patches. During the current equivalent calculation, it is essential that the source point and the field point are not obstructed by other patches. Otherwise, the current between the two must be ignored. In the IPO method, the occlusion judgment needs to be performed between any two patches. Due to the high computational complexity $ O\left( {{N^3}} \right) $ of the strict judgment process, the single occlusion judgment is often employed in the simulation of the IPO method. This involves judging the relationship between the incident direction $ {{\hat {\bf{k}}}_i} $ and the normal vector $ {\hat {\bf{n}}} $ of the patch. If $ {\hat {\bf{n}}} \cdot {{\hat {\bf{k}}}_i} < 0 $ is satisfied, it is determined to be a lit region. Otherwise, when $ {\hat {\bf{n}}} \cdot {{\hat {\bf{k}}}_i} > 0 $, it is determined to be a shadow region. However, the mutual occlusion between patches in complex target meshes can invalidate the results of this method to be invalid. Therefore, considering multiple occlusions could be beneficial to improving the iteration accuracy of the IPO method.

A strict shading method avoids false current equivalence, and on the other hand, it prevents additional current iterations as illustrated in Fig. 4. In short, a fast and rigorous shadow judgment method can enhance the scattering accuracy and tend to converge with fewer iterations, thereby reducing the simulation cost. Based on whether the ray and the patch number $ {\varDelta _m} $ are obstructed by other patches, determined by intersecting with the patches in the bounding box, if it is judged that the nearest patch number is the index $ m $, it is considered illuminated; otherwise, it is judged to be blocked or in shadow. Any point in the triangle $ {\varDelta _m} $ can be expressed as the position vector $ \overrightarrow {{\mathbf{AP}}} = u\overrightarrow {{\mathbf{AB}}} + v\overrightarrow {{\mathbf{AC}}} $, where the parameters $ u $ and $ v $ satisfy:

$ u = \frac{{\overrightarrow {{\mathbf{AQ}}} \cdot (\widehat {{\mathbf{QP}}} \times \overrightarrow {{\mathbf{AC}}} )}}{{\overrightarrow {{\mathbf{AB}}} \cdot (\widehat {{\mathbf{QP}}} \times \overrightarrow {{\mathbf{AC}}} )}}{\rm{,}}\; \tag{9a}$ ()

$v = \frac{{\overrightarrow {{\mathbf{AQ}}} \cdot (\widehat {{\mathbf{QP}}} \times \overrightarrow {{\mathbf{AB}}} )}}{{\overrightarrow {{\mathbf{AC}}} \cdot (\widehat {{\mathbf{QP}}} \times \overrightarrow {{\mathbf{AB}}} )}} \tag{9b}$ ()

where the $ {\mathbf{Q}} $ point originates from the triangle where the field point ${{\bf{r}}}$ is situated, with its midpoint typically selected. $ {\mathbf{A}} $, $ {{\bf{B}}} $, and $ {{\bf{C}}} $ represent the three endpoints of the triangle where the source point is located. Each position vector composed of these points can be observed in (9). If the point ${\bf{P}}(x{\rm{,}}\;y{\rm{,}}\;z)$ is inside the triangle, then $ 0 < u < 1 $, $ 0 < v < 1 $, and $ u + v < 1 $. After the ray intersects inside the triangle, the intersection distance between the ray and the triangle is

$ |{\mathbf{QP}}| = \left| {\frac{{\overrightarrow {{\mathbf{AC}}} \cdot (\overrightarrow {{\mathbf{AO}}} \times \overrightarrow {{\mathbf{AB}}} )}}{{\overrightarrow {{\mathbf{AB}}} \cdot (\widehat {{\mathbf{OP}}} \times \overrightarrow {{\mathbf{AC}}} )}}} \right|{\rm{.}}{\text{ }} $ (10)

In summary, an efficient shadow function judgment method for targets with multiple scattering contributions has been discussed. The high efficiency of the IPO-BVH method is verified with four numerical simulations in Section 5.

In this section, four numerical simulations are used to demonstrate the effectiveness of the proposed IPO-BVH algorithm. The case considered involves PEC in a far-field scenario. The simulation results are compared and verified against the outcomes obtained using the multi-level fast multipole algorithm in the commercial software FEldberechnung bei Körpern mit beliebiger Oberfläche (FEKO). Our proposed IPO-BVH method is showcased in numerical results. All simulations were performed in parallel on Intel Xeon Gold 6230R CPU with 64 GB RAM, 26 cores, and 52 threads.

An A380 aircraft inlet model was employed to verify the efficiency of the algorithm in the cavity target. When the target is of the same size, it is discretized into different numbers of patches. The CPU time consumption of the proposed IPO-BVH method and the initial method IPO-Brute in shadow judgment is depicted in Fig. 5. The number of patches $ M $ at the sampling points is distributed from 1374 to 108558. Each set of simulations selects the same one thousand patches to test the simulation time.

By comparing the traversal method of the original method with the method using the BVH structure, the results in Fig. 5 show that the proposed algorithm increases with the number of patches. The simulation time cost can be reduced by two orders of magnitude. The proposed method reduces the algorithm complexity of the shadow judgment process from $ O\left( {{N^3}} \right) $ to $ O\left( {{N^2}\;{\mathrm{log}}\; N} \right) $. When the number of patches $ M = 15710 $, the acceleration ratio of the proposed method is approximately 225. When $ M = 62572 $, the acceleration ratio rises to about 360. Since the method used avoids a large number of judgments through a tree structure, it exhibits excellent performance in the high-demand shadow occlusion judgment process.

The COBRA cavity model is selected for the accuracy verification of the proposed method. The radar cross section (RCS) patterns of higher-order finite element methods (FEM) for the COBRA model are demonstrated in Ref. [29]. The curing cavity was provided by European Aeronautic Defence and Space Company (EADS) Aerospatiale Matra Missiles for Workshop EM-JINA 98. Consider a plane wave with a simulation frequency of $ 17.5{\text{ }}{{\mathrm{GHz}}} $ to sweep the cavity aperture. The maximum number of iterations is set to 30 to determine the number of iterations on the inner wall of the cavity. The convergence threshold for each incident angle is set to 0.005. The lengths of the aperture in the ${x}$ and ${y}$ directions are $110{\text{ }}{{\mathrm{mm}}}$ and $84{\text{ }}{{\mathrm{mm}}}$, respectively. The length of the cavity in the $ {y} $ direction is $ 437{\text{ }}{{\mathrm{mm}}} $, and the length in the $ {z} $ direction is $ 128{\text{ }}{{\mathrm{mm}}} $. The cavity target is discretized into 13290 flat triangular patches, ensuring that each square contains at least $ 9 $ facets$ /{\lambda ^2} $. At this time, the incident pitch angle $ {\theta _{i}} = $ 90° and the azimuth angle $ {\varphi _{i}} = $ [90°, 150°]. A monostatic simulation was performed for horizontal polarization (H-H) and vertical polarization (V-V). The simulation results are illustrated in Fig. 6.

The MLFMA method is used as a standard solution to verify the results. The mean errors from the standard solution under V-V polarization and H-H polarization are $ 1.64{\text{ }}{{\mathrm{dB}}} $ and $ 1.33{\text{ }}{{\mathrm{dB}}} $, respectively. The equivalence principle at the deviation from the aperture plane will fail and lead to a decrease in the scattering accuracy of IPO [30]. Since the incident angle is close to the grazing aperture, the error from the standard solution becomes larger when the azimuth angle approaches $ {\varphi _{i}} = $ 150°. The stable distribution of the induced current under multiple incident azimuth angles $ {\varphi _{i}} = $ 90° and $ {\varphi _{i}} = $ 130° is shown in Fig. 7. In order to clearly display the current distribution inside the cavity, the current distribution on the virtual aperture surface is moved downward for illustration. The incident azimuth angle is $ {\varphi _{i}} = $ 90°, and the current is mainly distributed on the top and bottom surfaces of the cavity. As the incident angle increases to $ {\varphi _{i}} = $ 130°, the current distribution is mainly concentrated on the side of the cavity. The distribution of the current is related to the bounce of the ray inside the cavity.

The surface of the SLICY model has multiple scattering structures, including dihedral angles, angular reflections, and secondary scattering structures composed of cylinders and bottom surfaces. The virtual aperture surface $ {S_{a}} $ is not set here, so there is no need to obtain the initial iteration current of the target surface through the equivalence principle from (4). The target model is $ 0.61{\text{ }}{{\mathrm{m}}} $ long and $ 0.49{\text{ }}{{\mathrm{m}}} $ wide. The prism is a 1/4 cylinder with a radius of $ 61{\text{ }}{{\mathrm{mm}}} $, on which cylinders with diameters of $ 0.132{\text{ }}{{\mathrm{m}}} $ and $ 0.203{\text{ }}{{\mathrm{m}}} $ are placed, respectively, and the side length is $ 66{\text{ }}{{\mathrm{mm}}} $. The model is discretized into 16610 flat triangles and ensure that the density of 9 facets/$ {\lambda ^2} $. A plane wave with the simulated frequency $ f = 10{\text{ }}{{\mathrm{GHz}}} $ is incident along the pitch angle $ {\theta _{i}} $ = [0°, 90°] and the azimuth angle $ {\varphi _{i}} $ = 90°. The V-V and H-H polarization cases are considered in a monostatic case. The convergence threshold for current iterations is $ \varepsilon = 0.001 $. The mean errors of the results in the RCS pattern compared with the MLFMA method under the two polarization cases are $ 1.42{\text{ }}{{\mathrm{dB}}} $ and $ 2.31{\text{ }}{{\mathrm{dB}}} $, respectively. The results in Fig. 8 can prove that considering the strict shadow function between patches in Section 3 can ensure the better accuracy of the IPO algorithm and the full-wave algorithm.

The current distribution of IPO-BVH and PO methods for the plane wave with incident angle $ {\theta _{i}} $ = 45° and $ {\varphi _{i}} $ = 90° is shown in Fig. 9. Since the PO method only considers the approximate current of the initial iteration, its induced current distribution only depends on the normal vector. Compared with the current distribution of the IPO-BVH method, the distribution on the electrically large smooth-size target surface is smoother in Fig. 9 (a). However, at incident angles with multiple scattering contributions, the accuracy of PO drops significantly from Fig. 8. However, under the strict judgment of the patch shadow relationship in (10) at this time, the proposed method can still maintain well consistency with the standard solution.

In order to testify the applicability of the algorithm in general scenarios, the ship-sea scenary is considered for verification. On the rough surface, where the length and width measure $8{\text{ }}{{\mathrm{m}}}$, there is a ship with dimensions of $7.4{\text{ }}{{\mathrm{m}}}$ in length, $0.8{\text{ }}{{\mathrm{m}}}$ in width, and $1.6{\text{ }}{{\mathrm{m}}}$ in height. The rough surface model is generated with the Monte Carlo method based on the Pierson and Moskowitz (PM) spectrum. The ship is in a draft state on the sea model. The scene is discretely generated into 19252 triangular patches according to $ 1/3\lambda $ wavelength. The V-V polarization at the incident plane wave azimuth angle $ {\varphi _{i}} $ = 0° and $ {\varphi _{i}} $ = 90°, respectively, and the pitch angle $ {\theta _{i}} $ = [0°, 70°] is considered. In the case of the monostatic station, and the threshold is set to 0.005. The simulation results are shown in Fig. 10, and the mean errors with the MLFMA results are $ 1.14{\text{ }}{{\mathrm{dB}}} $ and $ 3.12{\text{ }}{{\mathrm{dB}}} $, respectively. When the incident azimuth angle is $ {\varphi _{i}} $ = 0°, as the pitch angle gradually increases, the number of iterations gradually increases from 2 to 4.

The current distribution of PO and IPO-BVH methods is shown in Fig. 11. The significant difference in amplitudes of induced current between the PO method and the IPO-BVH method arises from the contribution of multiple scattering. Due to the angular inverse structure formed between the ship and the sea surface, the IPO-BVH method can maintain the high accuracy after considering the strict shadow function between patches in Section 3. This can explain that as the incident angle gradually tends to a low grazing angle, the iteration between the ship and the sea surface gradually changes to the iteration between ships. On the other hand, this method is suitable for judging the shadow function between patches in complex scenes. However, judging only by normal vectors will lead to the incorrect current iteration between patches [19]. Therefore, the proposed method is beneficial to the convergence process of the current iteration.

In this paper, BVH is proposed to accelerate the shadow judgment process between dense patches in the IPO methods. The BVH tree of the model is constructed by constructing a binary tree, thus avoiding the complexity of $ O\left( {{N^3}} \right) $ in the traditional ray intersection method and reducing it to $ O\left( {{N^2}\;{\mathrm{log}}\; N} \right) $. Numerical simulations prove that the proposed method can improve the simulation accuracy and achieve high efficiency in the obtaining simulation results. Numerical results show that the proposed method can reach the convergence threshold. In the future, we will also apply the proposed method to the IPO method based on high-order patch unit discretization and need to consider the problem of the iteration non-convergence in some scenarios.

The authors declare no conflicts of interest.