# The Efficient Method for Calculating the Physical Optics Scattered Fields from the Concave Surfaces

Authors: Yu Mao Wu; Hanzhang Zhou; Ya-Qiu Jin; Jun Hu; Haijing Zhou; Weijie Wang
Journal: International Journal of Antennas and Propagation
 (2020)

---

## Abstract

In this work, the numerical steepest descent path (NSDP) method is proposed to compute the highly oscillatory physical optics (PO) scattered fields from the concave surfaces, including both the monostatic and the bistatic cases. Quadratic variations are adopted to approximate the integrands of the PO type integral into the canonical form. Then, on involving the NSDP method, we deform the integration paths of the integrals into several NSDPs on the complex plain, through which the highly oscillatory integrands are converted to exponentially decay integrands. The RCS results of the PO scattered field are calculated and are compared with the high frequency asymptotic (HFA) method and the brute force (BF) method. The results demonstrate that the proposed NSDP method for calculating PO scattered fields from concave surfaces is frequency-independent and error-controllable. Numerical examples are provided to verify the efficiencies of the NSDP method.


---

## Body

## 1. Introduction

In electromagnetics (EM) community, as the product of the wavenumberk and the electrical size of the considered object d are large enough, the calculation of the scattered EM field belongs to high frequency scattering problems. Computing the high frequency scattered fields from electrically large perfect conducting scatterer remains significant and challenging [1–4]. In 1913, Macdonald introduced the physical optics (PO) approximation [5], which is used as an effective method for computing EM scattering problems [4, 6–8]. By this method, local tangent plane approximation is adopted to approximate the current in lit parts, and current in shadow parts is taken as zero. The PO approximation [5, 6, 8–14] is of high efficiency when calculating scattered fields from smooth and electrically large scatterer. Mathematically, the amplitude term of PO integrand is slowly varying. However, the exponential of the phase function term is highly oscillatory as the frequency k increases. In this sense, the computational effort of a direct numerical method [2, 11] for calculating the PO type integral is extremely high.Some studies assume that the amplitude function and phase function terms vary linearly. Then, the physical optics integral is analytically simplified into several line integrals [15–17]. In this situation, the closed form formulas can be derived from the flat polygon patches. However, the stationary phase point (SPP), the resonance points (RP), and the vertex point (VP) contributions of the PO integral can not be captured when the phase terms of the PO type integrand are approximated into linear terms. Therefore, the traditional high frequency asymptotic (HFA) method [6, 18–21] was proposed, approximating the phase terms of the integral into quadratic forms. The HFA method could approximate the PO integrand into several dominating terms [18], corresponding to contributions from the SPP, the RP, and the VP. With the aid of HFA method, the PO-based EM fields could be calculated with frequency-independent computational workload.The numerical steepest descent path (NSDP) approach [22–29] is proposed as an effective way to calculate the highly oscillatory PO type integral. With the aid of the NSDP method, the highly oscillatory PO type integrand is converted into several smooth ones on the complex plane by the contour deformation technique. Importantly, the contributions from the SPP, the RP, and the VP were captured through the NSDP method successfully. In former works [26, 27], the PO type integral from hyperbolic and quadratic convex surfaces has been studied. Our recent study [29] summarizes the promising progress on calculating high frequency scattered fields with the aid of the NSDP method. However, basic theories on NSDP approach have not been completed yet. Specifically, the NSDP method needs to be extended to the realistic concave surfaces.The contributions of this work are that, first, physical optics scattered fields from concave surfaces are considered. In this case, the phase function term of the PO integrand takes the quadratic concave behavior, which differs from works in [26, 27]. Therefore, the derivations of calculating the PO type integral are different, the NSDPs are correspondingly changed, the stationary phase points and boundary resonance points are all different, and the integral paths are also changed. Furthermore, by adopting the NSDP method, the PO type integral from concave surfaces could be calculated in frequency-independent workload and error-controllable accuracy. Compared with the HFA method, the accuracy could be improved by around 2 digits. Finally, former work in [28] only calculates the PO type integral from limited concave triangular patches. In this work, however, we extend the NSDP method to the realistic concave surfaces. Together with works in [26, 27], the NSDP method could be employed to treat PO integral with all types of phase behaviors and be applied to realistic scatterer.This communication is organized as follows. First, the PO integral under the PO approximation is derived, and the considered concave surface is discretized into several triangular patches. The Lagrange interpolation polynomial approximation and the affine transformation are used to approximate the amplitude function and phase function of the PO type integral. Next, the double integral is transformed to some highly oscillatory line integrals after integration by part. Then, we propose the NSDP method to obtain corresponding NSDPs and convert highly oscillatory PO type integrand into several smooth ones. Finally, numerical examples of calculating PO scattered fields from concave surfaces are provided, and comparisons with the brute force method and HFA method are given.


## 2. PO Surface Integral

### 2.1. Physical Optics Approximation

In this work, the size of considered concave surface is denoted asL=diam∂Ω, the distance between the observation point r and the origin is r=r^2, and wavelength is λ. We consider the far-zone of the scatterer with r≫2L2/λ and adopt the PO approximation to the accurate electric equivalent current. Then, according to the Stratton–Chu integral equation formulation [3], for the plane incident wave, the PO scattered fields can be written as(1)Esr′,r≈∫∂Ωsbir′,reikvbir′,rdSr′,(2)sbir′,r=−ikeikr2πrr^×r^×n^r′×r^i×E0i,(3)vbir′,r=r′^−r^⋅r′.Here,k is the wavenumber outside ∂Ω, E0i is the incident electric polarization wave vector, r is the observation point, and r′ is the surface point on ∂Ω, as shown in Figure 1(a).Figure 1
(a) The 3D illustration of the quadratic concave triangular patch. (b) The projection of the concave triangular patch on thex-y plane.
(a)(b)Equation (1) is the expression of the PO type integral of bistatic scattered electric field. The vector function sbir′ in (2) is the vector amplitude function and vbir′ in (3) is the phase function.For the monostatic scattered field, withr^=−r^i, (1) can be expressed as(4)Esr′,r≈∫∂Ωsmor′,reikvmor′,rdSr′,with(5)smor′,r=−ikeikr2πrr^i⋅n^r′,vmor′,r=2r^i⋅r′.Therefore,Esr′,r under the PO approximation for both bistatic and monostatic cases takes the general form [25].(6)I˜=∫∂Ωsr′,reikvr′,rdSr′.The integrand contains a slowly varying amplitude term and the exponential of the phase term. Whenk increases, the integrand is highly oscillatory due to the eikvr′,r term, as shown in Figures 2(a)–2(d). In this sense, the computational cost of a direct numerical method grows dramatically with a larger k.Figure 2
(a, b) The real and imag part of a PO type integrand functionFx=6−x2eikx2+x+12 when k=1. (c, d) The real and imag part of Fx when k=100. (e, f) The real and imag part of Fϕ0p when k=100, which is derived by defining Fx on contour path ϕ0p.
(a)(b)(c)(d)(e)(f)


### 2.2. Extending the Triangular Patch into the Realistic Quadratic Surface

In previous work [28], we preliminarily calculated the PO type integral from quadratic concave surfaces. However, the method had clear limitation; that is, it could only be applied for triangular patches and was not able to calculate the realistic quadratic surfaces. On the other hand, the actual geometric modeling is generally carried out by quadratic surface subdivision. In this sense, the previous work has been flawed in engineering applications, and it is necessary to develop numerical method for the realistic quadratic surfaces.We first project the quadratic surfaceΩ into the xy coordinate system, which is Ωxy. Then, Ωxy is divided into several triangular patches. In this way, the PO surface integral (7) can be expressed as(7)I=∫∂Ωsx,y,zeikvx,y,zdxdydz=∫∂Ωxys˜x,yeikv˜x,y1+fxx,y2+fyx,y2dydx=∑n=1M∫Δnd˜nx,yeikv˜nx,ydydx.Here,z=fx,y is the expression of the quadratic surface Ω. The second-order polynomials d˜nx,y and v˜nx,y in (7) are the approximated amplitude and phase functions, which are obtained by the Lagrange interpolation polynomial approximation [30] on the triangular patches Δn,n=1,2,…,M, respectively. Formulas of d˜nx,y and v˜nx,y are(8)v˜nx,y=β˜n,1+β˜n,2x+β˜n,3y+β˜n,4x2+β˜n,5y2+β˜n,6xy,d˜nx,y=α˜n,1+α˜n,2x+α˜n,3y+α˜n,4x2+α˜n,5y2+α˜n,6xy,with coefficients β˜n,m∈ℂ,α˜n,m∈ℝ,m=1,2,…,6,n=1,2,…,M.In addition, in order to simplify the derivation, affine transformation is used to simplify the phase functionv˜nx,y into its canonical form:(9)I=∑n=1M∫Δn′p˜nx′,y′eik±x′2±y′2dy′dx′.


## 2.1. Physical Optics Approximation

In this work, the size of considered concave surface is denoted asL=diam∂Ω, the distance between the observation point r and the origin is r=r^2, and wavelength is λ. We consider the far-zone of the scatterer with r≫2L2/λ and adopt the PO approximation to the accurate electric equivalent current. Then, according to the Stratton–Chu integral equation formulation [3], for the plane incident wave, the PO scattered fields can be written as(1)Esr′,r≈∫∂Ωsbir′,reikvbir′,rdSr′,(2)sbir′,r=−ikeikr2πrr^×r^×n^r′×r^i×E0i,(3)vbir′,r=r′^−r^⋅r′.Here,k is the wavenumber outside ∂Ω, E0i is the incident electric polarization wave vector, r is the observation point, and r′ is the surface point on ∂Ω, as shown in Figure 1(a).Figure 1
(a) The 3D illustration of the quadratic concave triangular patch. (b) The projection of the concave triangular patch on thex-y plane.
(a)(b)Equation (1) is the expression of the PO type integral of bistatic scattered electric field. The vector function sbir′ in (2) is the vector amplitude function and vbir′ in (3) is the phase function.For the monostatic scattered field, withr^=−r^i, (1) can be expressed as(4)Esr′,r≈∫∂Ωsmor′,reikvmor′,rdSr′,with(5)smor′,r=−ikeikr2πrr^i⋅n^r′,vmor′,r=2r^i⋅r′.Therefore,Esr′,r under the PO approximation for both bistatic and monostatic cases takes the general form [25].(6)I˜=∫∂Ωsr′,reikvr′,rdSr′.The integrand contains a slowly varying amplitude term and the exponential of the phase term. Whenk increases, the integrand is highly oscillatory due to the eikvr′,r term, as shown in Figures 2(a)–2(d). In this sense, the computational cost of a direct numerical method grows dramatically with a larger k.Figure 2
(a, b) The real and imag part of a PO type integrand functionFx=6−x2eikx2+x+12 when k=1. (c, d) The real and imag part of Fx when k=100. (e, f) The real and imag part of Fϕ0p when k=100, which is derived by defining Fx on contour path ϕ0p.
(a)(b)(c)(d)(e)(f)


## 2.2. Extending the Triangular Patch into the Realistic Quadratic Surface

In previous work [28], we preliminarily calculated the PO type integral from quadratic concave surfaces. However, the method had clear limitation; that is, it could only be applied for triangular patches and was not able to calculate the realistic quadratic surfaces. On the other hand, the actual geometric modeling is generally carried out by quadratic surface subdivision. In this sense, the previous work has been flawed in engineering applications, and it is necessary to develop numerical method for the realistic quadratic surfaces.We first project the quadratic surfaceΩ into the xy coordinate system, which is Ωxy. Then, Ωxy is divided into several triangular patches. In this way, the PO surface integral (7) can be expressed as(7)I=∫∂Ωsx,y,zeikvx,y,zdxdydz=∫∂Ωxys˜x,yeikv˜x,y1+fxx,y2+fyx,y2dydx=∑n=1M∫Δnd˜nx,yeikv˜nx,ydydx.Here,z=fx,y is the expression of the quadratic surface Ω. The second-order polynomials d˜nx,y and v˜nx,y in (7) are the approximated amplitude and phase functions, which are obtained by the Lagrange interpolation polynomial approximation [30] on the triangular patches Δn,n=1,2,…,M, respectively. Formulas of d˜nx,y and v˜nx,y are(8)v˜nx,y=β˜n,1+β˜n,2x+β˜n,3y+β˜n,4x2+β˜n,5y2+β˜n,6xy,d˜nx,y=α˜n,1+α˜n,2x+α˜n,3y+α˜n,4x2+α˜n,5y2+α˜n,6xy,with coefficients β˜n,m∈ℂ,α˜n,m∈ℝ,m=1,2,…,6,n=1,2,…,M.In addition, in order to simplify the derivation, affine transformation is used to simplify the phase functionv˜nx,y into its canonical form:(9)I=∑n=1M∫Δn′p˜nx′,y′eik±x′2±y′2dy′dx′.


## 3. Calculating the PO Scattered Fields With the Concave Phase Function Term

### 3.1. The Concave Phase Function Case

From (9), it can be concluded that there are four different cases of the phase function due to different combination of phase term, which are x2+y2, x2−y2, −x2+y2, and −x2−y2, respectively. Calculating PO surface scattered field with phase function of the first three cases has been studied by our previous works [26, 27]. In this work, we study the concave surfaces case, and a canonical PO type integral can be derived from (9):(10)In=∫Δnpx,ye−ikx2+y2dydx.With different phase function term of the PO integrand, derivations of calculating the integral shall be changed. Most importantly, the numerical steepest decent paths shall all be changed, and the stationary phase points and boundary resonance points are all different. As a result, the spacial relationships between stationary phase points and end points and numerical steepest decent paths and stokes line are all correspondingly changed. Therefore, the integral paths after adopting contour deformation technique are also changed. The NSDP method for the concave surface is proposed to compute the canonical PO type integralIn in a frequency-independent workload.


### 3.2. Reduction of the Double Integral into Line Integral

The amplitude term of (10), px,y, has the expression(11)px,y=α1+α2x+α2y+α4x2+α5y2+α6xy.Then, after integration by part, we obtain(12)I=∫L1L2∫ax+b0px,ye−ikx2+y2dydx=∫L1L2∫ax+b0α1+α2x+α4x2e−ikx2+y2dydx+∫L1L2∫ax+b0α2y+α6xye−ikx2+y2dydx+∫L1L2∫ax+b0α5y2e−ikx2+y2dydx=∫L1L2F1xe−ikx2dx+∫L1L2F2xe−ikx2dx+∫L1L2F3xe−ikx2dx.Then(13)Fx=F1x+F2x+F3x=−j1xerfcikax+b−j2a,bxe−ikax+b2+j20,0x+j1x,with(14)j1x=−π2ikα1+α2x+α4x2+α52ik,j2a,bx=α3+α6x+α5ax+b2ki,j20,0x=α3+α6x2ki.In (13), there are two terms defined on the edge y=ax+b and another two terms defined on the edge y=0, and the term defined on y=0 is easy to calculate. Therefore, Fx is divided into J2a,bx and J20,0x. They can be expressed as(15)J20,0x=j1x+j20,0x,J2a,bx=j1xerfcikax+b+j2a,bxe−ikax+b2.Hence, we can rewriteI as(16)I=∫L1L2−J2a,bxe−ikx2dx+∫L1L2J20,0xe−ikx2dx=−I2+I1.In (15), both I1 and I2 are highly oscillatory integrals, but I1 has the form of complementary error function and can be easily calculated with the aid of complementary error function. I2, however, can not be calculated analytically and the computational cost increases dramatically with the larger k. In the following, numerical steepest descent path method is illustrated to calculate I2.


### 3.3. Calculating the Highly Oscillatory IntegralI2

The contour deformation path technique is adopted to calculateI2. After substituting y=ax+b, we could obtain the phase function term of the integrand for I2 as(17)gx=−x2−ax+b2.The proposed steepest descent path should satisfy the following three conditions,p∈0,∞:(a)
φm0=Lm, which means that the NSDPs start from the end points Lm, m=1,2(b)
Regφmp=Regφm0≡C, where C is a constant(c)
Imgφmp=pIn this sense, the exponential term of the phase function can be derived as(18)eikgx=eikRegx+iImgx=e−kImgx+ikRegx=e−kp+ikC.It can be seen thateikC term is a certain imaginary number and ekp is a real term that decreases exponentially when p goes large. Therefore, the exponential of phase function is clearly no longer oscillatory. The contrast between Figures 2(c)–2(d) and Figures 2(e)–2(f) illustrates the effect of our NSDP method.From the definition of NSDPs, we could obtain steepest descent path for end points asφmp, m=1,2:(19)φmp=signLm′1+a2Lm′2−ip+xs,p∈0,∞,where(20)Lm′=a2+1Lm−xs,xs=−aba2+1,where xs is the stationary phase point (SPP), which satisfies g′xs=0. The phase behavior of the SPP differs from the end points. In this sense, we define a contour path function for SPP as x=φ0p,p∈−∞,∞:(21)gφ0p=gxs+ip2.We could also obtain NSDP forxs as follows:(22)φ0p=e−iπ/4p1+a2+xs,p∈−∞,∞.By adopting the NSDP method,I2 can be written as(23)I2=∫φNSDPpj1xerfcikax+b+j2a,bxe−ikax+b2⋅e−ikx2dx.At this point, the behavior of complementary error function should be considered.erfcikz has the following feature:(24)erfcikz=ζ1ze−ikz2,Rez−Imz≤0,2+ζ2ze−ikz2,Rez−Imz>0,where ζ1z and ζ2z are two slowly varying functions.The above different behavior oferfcikz on two different regions is called Stokes phenomenon. The line separating these two regions is called Stokes line, which can be expressed as(25)Rez−Imz=0.After substitutingz=ax+b into (25), the Stokes line takes the following expression:(26)Rex=Imx−ab.Considering the features of erfcikz,I2 in (23) can be further written as(27)I2=I2non−analytic+I2analytic,where(28)I2non−analytic=∫φNSDPj1xζx+j2a,bxeikgxdx,I2analytic=∫φ∗2j1xe−ikx2dx,ζx=ζ1ax+b,aRex−Imx≤−b,ζ2ax+b,aRex−Imx>−b.For integralI2analytic, the original function of its integrand is(29)Kx=π2ikα1−π4k2α4−π4k2α5erfcikx+π2ikikα2+πx2ikikα4e−ikx2,and I2analytic is from the constant “2” in (24), and its integral path φ∗ is determined by the spacial relationship between Stokes line and NSDPs. Specifically, different cases of integral paths arise from relative locations of end points, SPP, and intersection point of Stokes line and x axis. Then I2analytic can be calculated.ForI2non−analytic, its specific expression comes from the relative values of end points L1,L2, and SPP xs, which can be obtained by the following equations:(30)I2non−analytic=I2φ1−I2φ2+I2φ0,L1<xs<L2,I2φ1−I2φ2,L2<xsorxs<L1.After substitutingx with the corresponding path function, we could obtain the expression for the terms in (30):(31)I2φ1=∫0∞e−kpe−ikL12+aL1+b2K1pdp,(32)I2φ2=∫0∞e−kpe−ikL22+aL2+b2K2pdp,(33)I2φ0=∫−∞∞e−kp2e−ikxs2+axs+b2K0pdp.Here,(34)K1p=j1φ1pζφ1p+j2a,bφ1pφ1′p,K2p=j1φ2pζφ2p+j2a,bφ2pφ2′p,K0p=j1φ0pζφ0p+j2a,bφ0pφ0′p.At this point,I2analytic and I2non−analytic are obtained, and I2 can be calculated efficiently.


## 3.1. The Concave Phase Function Case

From (9), it can be concluded that there are four different cases of the phase function due to different combination of phase term, which are x2+y2, x2−y2, −x2+y2, and −x2−y2, respectively. Calculating PO surface scattered field with phase function of the first three cases has been studied by our previous works [26, 27]. In this work, we study the concave surfaces case, and a canonical PO type integral can be derived from (9):(10)In=∫Δnpx,ye−ikx2+y2dydx.With different phase function term of the PO integrand, derivations of calculating the integral shall be changed. Most importantly, the numerical steepest decent paths shall all be changed, and the stationary phase points and boundary resonance points are all different. As a result, the spacial relationships between stationary phase points and end points and numerical steepest decent paths and stokes line are all correspondingly changed. Therefore, the integral paths after adopting contour deformation technique are also changed. The NSDP method for the concave surface is proposed to compute the canonical PO type integralIn in a frequency-independent workload.


## 3.2. Reduction of the Double Integral into Line Integral

The amplitude term of (10), px,y, has the expression(11)px,y=α1+α2x+α2y+α4x2+α5y2+α6xy.Then, after integration by part, we obtain(12)I=∫L1L2∫ax+b0px,ye−ikx2+y2dydx=∫L1L2∫ax+b0α1+α2x+α4x2e−ikx2+y2dydx+∫L1L2∫ax+b0α2y+α6xye−ikx2+y2dydx+∫L1L2∫ax+b0α5y2e−ikx2+y2dydx=∫L1L2F1xe−ikx2dx+∫L1L2F2xe−ikx2dx+∫L1L2F3xe−ikx2dx.Then(13)Fx=F1x+F2x+F3x=−j1xerfcikax+b−j2a,bxe−ikax+b2+j20,0x+j1x,with(14)j1x=−π2ikα1+α2x+α4x2+α52ik,j2a,bx=α3+α6x+α5ax+b2ki,j20,0x=α3+α6x2ki.In (13), there are two terms defined on the edge y=ax+b and another two terms defined on the edge y=0, and the term defined on y=0 is easy to calculate. Therefore, Fx is divided into J2a,bx and J20,0x. They can be expressed as(15)J20,0x=j1x+j20,0x,J2a,bx=j1xerfcikax+b+j2a,bxe−ikax+b2.Hence, we can rewriteI as(16)I=∫L1L2−J2a,bxe−ikx2dx+∫L1L2J20,0xe−ikx2dx=−I2+I1.In (15), both I1 and I2 are highly oscillatory integrals, but I1 has the form of complementary error function and can be easily calculated with the aid of complementary error function. I2, however, can not be calculated analytically and the computational cost increases dramatically with the larger k. In the following, numerical steepest descent path method is illustrated to calculate I2.


## 3.3. Calculating the Highly Oscillatory IntegralI2

The contour deformation path technique is adopted to calculateI2. After substituting y=ax+b, we could obtain the phase function term of the integrand for I2 as(17)gx=−x2−ax+b2.The proposed steepest descent path should satisfy the following three conditions,p∈0,∞:(a)
φm0=Lm, which means that the NSDPs start from the end points Lm, m=1,2(b)
Regφmp=Regφm0≡C, where C is a constant(c)
Imgφmp=pIn this sense, the exponential term of the phase function can be derived as(18)eikgx=eikRegx+iImgx=e−kImgx+ikRegx=e−kp+ikC.It can be seen thateikC term is a certain imaginary number and ekp is a real term that decreases exponentially when p goes large. Therefore, the exponential of phase function is clearly no longer oscillatory. The contrast between Figures 2(c)–2(d) and Figures 2(e)–2(f) illustrates the effect of our NSDP method.From the definition of NSDPs, we could obtain steepest descent path for end points asφmp, m=1,2:(19)φmp=signLm′1+a2Lm′2−ip+xs,p∈0,∞,where(20)Lm′=a2+1Lm−xs,xs=−aba2+1,where xs is the stationary phase point (SPP), which satisfies g′xs=0. The phase behavior of the SPP differs from the end points. In this sense, we define a contour path function for SPP as x=φ0p,p∈−∞,∞:(21)gφ0p=gxs+ip2.We could also obtain NSDP forxs as follows:(22)φ0p=e−iπ/4p1+a2+xs,p∈−∞,∞.By adopting the NSDP method,I2 can be written as(23)I2=∫φNSDPpj1xerfcikax+b+j2a,bxe−ikax+b2⋅e−ikx2dx.At this point, the behavior of complementary error function should be considered.erfcikz has the following feature:(24)erfcikz=ζ1ze−ikz2,Rez−Imz≤0,2+ζ2ze−ikz2,Rez−Imz>0,where ζ1z and ζ2z are two slowly varying functions.The above different behavior oferfcikz on two different regions is called Stokes phenomenon. The line separating these two regions is called Stokes line, which can be expressed as(25)Rez−Imz=0.After substitutingz=ax+b into (25), the Stokes line takes the following expression:(26)Rex=Imx−ab.Considering the features of erfcikz,I2 in (23) can be further written as(27)I2=I2non−analytic+I2analytic,where(28)I2non−analytic=∫φNSDPj1xζx+j2a,bxeikgxdx,I2analytic=∫φ∗2j1xe−ikx2dx,ζx=ζ1ax+b,aRex−Imx≤−b,ζ2ax+b,aRex−Imx>−b.For integralI2analytic, the original function of its integrand is(29)Kx=π2ikα1−π4k2α4−π4k2α5erfcikx+π2ikikα2+πx2ikikα4e−ikx2,and I2analytic is from the constant “2” in (24), and its integral path φ∗ is determined by the spacial relationship between Stokes line and NSDPs. Specifically, different cases of integral paths arise from relative locations of end points, SPP, and intersection point of Stokes line and x axis. Then I2analytic can be calculated.ForI2non−analytic, its specific expression comes from the relative values of end points L1,L2, and SPP xs, which can be obtained by the following equations:(30)I2non−analytic=I2φ1−I2φ2+I2φ0,L1<xs<L2,I2φ1−I2φ2,L2<xsorxs<L1.After substitutingx with the corresponding path function, we could obtain the expression for the terms in (30):(31)I2φ1=∫0∞e−kpe−ikL12+aL1+b2K1pdp,(32)I2φ2=∫0∞e−kpe−ikL22+aL2+b2K2pdp,(33)I2φ0=∫−∞∞e−kp2e−ikxs2+axs+b2K0pdp.Here,(34)K1p=j1φ1pζφ1p+j2a,bφ1pφ1′p,K2p=j1φ2pζφ2p+j2a,bφ2pφ2′p,K0p=j1φ0pζφ0p+j2a,bφ0pφ0′p.At this point,I2analytic and I2non−analytic are obtained, and I2 can be calculated efficiently.


## 4. An Example of Calculating PO Surface Integral by Adopting the NSDP Method

We consider a quadratic concave triangular patchΔV1′V2′V3′, as shown in Figure 1(a). A flat triangular patch with three vertex points V1=V1x,V1y, V2=V2x,V2y, V3=V3x,V3y could be derived by projecting ΔV1′V2′V3′ onto the x-y plane, as shown in Figure 1(b). And the three edges of ΔV1V2V3 are governed by y=amx+bm; the phase functions on these edges are gmx=−x2−amx+bm2, m=1,2,3. We define V1V2, V2V3, and V3V1 as edge 1, edge 2, and edge 3, respectively. As I1 can be calculated analytically, in the following, we will discuss I2 on three different edges of the triangular patch, respectively. First, taking the integral on edge 1 V1V2 as an example, its steepest descent paths and Stokes line are shown in Figure 3(a).Figure 3
Steepest descent paths for the integrand ofI2 defined on edge 1 (a), edge 2 (b), and edge 3 (c).
(a)(b)(c)For edge 1, which isV1V2, V1x, V2x, and y=a1x+b1 correspond to L1, L2, and y=ax+b in the last section, respectively. Then, I in (16) takes the following expression:(35)Iedge1=I2,a1,b1non−analytic+I2,a1,b1analytic−I1,a1,b1.As in Figure3(a), V1x<V2x<xs; according to (30), I2,a1,b1non−analytic can be written as(36)I2,a1,b1non−analytic=I2,a1,b1V1x−I2,a1,b1V2x,where I2,a1,b1V1x and I2,a1,b1V2x have the similar forms as I2φ1 and I2φ2 given in (31).As forI2,a1,b1analytic, B=B1x,B1y is the intersection point of Stokes line and steepest descent paths. Considering the integral path φ∗ is in the domain satisfying aRex−Imx>−b, which is the domain above Stokes line, we can obtain I2,a1,b1analytic as follows:(37)I2,a1,b1analytic=∫V1xB1x2j1xe−ikx2d=KB1x−KV1x,where KB1x and KV1x have the similar forms as Kx given in (30). At this point, I2 for edge 1 is obtained.(38)I2,edge1=I2,a1,b1V1x−I2,a1,b1V2x+KB1x−KV1x.Similarly, for edge 2,V2V3 (L1=V2x, L2=V3x, and y=a2x+b2), we have(39)Iedge2=I2,edge2−I1,a2,b2,I2,edge2=I2,a2,b2V2x+I2,a2,b2xs2−I2,a2,b2V3x+KV3x−KA2x.To be noticed, as shown in Figure3(b), V1x<xs<V2x. Therefore, according to (29), I2,a2,b2non−analytic contains three terms, which are I2,a2,b2V2x, I2,a2,b2xs2, and I2,a2,b2V3x, respectively.And for edge 3,V1V3 (L1=V1x, L2=V3x, and y=a3x+b3), we have(40)Iedge3=I2,edge3−I1,a3,b3.Here,(41)I2,edge3=I2,a3,b3V1x+I2,a3,b3xs3−I2,a3,b3V3x+KC3x−KV1x+KV3x−KA3x.Based on (42)–(44), we could derive the closed form formulas for I defined on quadratic concave triangular patches:(42)I=Iedge1+Iedge2−Iedge3=I2,a1,b1V1x−I2,a3,b3V1x+I2,a2,b2V2x−I2,a1,b1V2x+I2,a3,b3V3x−I2,a2,b2V3x+I2,a2,b2xs2−I2,a3,b3xs3+KB1x−KA2x+KA3x−KC3x,where I1 is canceled during the NSDP technique process. We can also separate the contributions of three critical points, which are the stationary phase point (SPP), the resonance points (RP), and the vertex points (VP):(43)IVP=I2,a1,b1V1x−I2,a3,b3V1x+I2,a2,b2V2x−I2,a1,b1V2x+I2,a3,b3V3x−I2,a2,b2V3x,(44)IRP=I2,a2,b2xs2−I2,a3,b3xs3,(45)ISPP=KB1x−KA2x+KA3x−KC3x.Hence, the PO surface integralI can be expressed by the high frequency critical point contributions as follows:(46)I=IVP+IRP+ISPP.


## 5. Numerical Results

### 5.1. Calculating the Assembled Concave Quadratic Triangular Patches

We consider a concave quadratic patch (see Figure4):(47)fx,y,z:z=1+0.6x2+xy+y2,with four vertex points: V1=0.8932,−13.3333, V2=22.7671,−16.6667, V3=9.8803,3.3333, and V4=−11.2201,16.6667. The incident wave direction ri=0.5,0.5,−2/2 with the amplitude Ei=−0.5,0.5,0 are considered when we treat a bistatic situation. The far field observation point lies along the direction r¯s=2/4,6/4,2/2T. The quadrilateral domain V1V2V3V4 in Figure 4(b) is decomposed into two triangular domains, ΔV1V2V3 and ΔV1V3V4, respectively.Figure 4
The original quadratic patch. (a) 3D illustration. (b) 2D illustration.
(a)(b)On invoking the affine transformation, we obtain the transformed new quadratic patch. The new vertex points for the transformed new patch areV˜1=−1.9331,−2.9957,V˜2=3.9531,−3.8366,V˜3=3.1122,1.2087,V˜4=−1.0923,4.5723.The concave quadratic patch can be calculated through two assembled triangular patches based on the NSDP method introduced above. The computational results and the corresponding consumed CPU time with increasing wavenumberk are also calculated by the high frequency asymptotic (HFA) method and the brute force (BF). Specifically, from Figure 5(a), it can be inferred that the computational results of the NSDP method agree well with those from BF method and the HFA method. Furthermore, Figure 5(b) demonstrates that the NSDP method is able to calculate the scattered fields from concave surfaces in a frequency-independent manner.Figure 5
Results and CPU time of the scattered fields with the increasing frequencyk. (a) Comparisons of the RCS values produced by the NSDP method, the BF method, and the HFA method. (b) Comparisons of the CPU time of calculating the PO scattered fields by the NSDP method, the BF method, and the HFA method.
(a)(b)


### 5.2. The Efficiency of the NSDP Method

We consider another concave PEC patch:(48)fx,y,z:z=1+0.06x2+xy+y2,to benchmark the efficiencies of our proposed NSDP method on calculating the PO scattered fields.The projection of the surface on thexy plane is a quadrilateral domain, as shown in Figure 6, with four vertex points: V1=9.0069,−13.3791, V2=22.7671,−16.6667, V3=9.8803,3.3333, and V4=−11.2201,16,6667. After the affine transformation, the quadrilateral domain in Figure 6(a) is transformed into another quadrilateral domain in Figure 6(b). We set the parameters as follows: the wave frequency k∈0,400 and the incident wave propagates along r^i=0.5,0.5,−2/2T direction. Besides the NSDP method, both the BF method and the HFA method are applied to calculate the PO scattered fields to illustrate the efficiency of the NSDP method.Figure 6
(a) The projection of the concave surface on thexy plane: V1V2V3V4. (b) The projection VV1VV2VV3VV4 after the affine transformation.
(a)(b)Through the BF method verification, we could compare the results generated by the NSDP method and the HFA method. The computational results of the PO scattered field and the relative error generated by the NSDP and the HFA methods are shown in Table1 and Table 2, respectively. Meanwhile, the CPU time consumed by these methods is given in Table 3. Figure 7(a) illustrates the RCS values of the PO scattered field Esr produced by both methods, which match well with results from the BF method. Moreover, Figure 7(b) indicates that the accuracy can be significantly improved by around two digits 102 through the NSDP method. Figure 7(c) demonstrates the CPU time consumed by both the NSDP method and the HFA method. Clearly, the computational workload for the NSDP method is relatively low and frequency-independent.Table 1
Computational results of the RCS values (dBsm unit) of PO scattered fields from the concave patch by using the BF method, the NSDP method, and the HFA method.
FrequencykBF-RCSNSDP-RCSHFA-RCS2024.748124.756724.72945026.012226.010825.986810026.102326.110726.233930027.465627.499927.437740024.917624.924524.9925Table 2
Relative error (using the results of the BF method as standard values) of calculating the PO scattered fields by using the NSDP method and the HFA method.
FrequencykNSDPHFABF203.4478×10−44.0559×10−2—505.1694×10−59.7626×10−4—1003.2366×10−45.0415×10−3—3008.8057×10−45.4585×10−3—4006.3527×10−55.0082×10−3—Table 3
Consumed CPU time (second unit) of calculating the PO scattered fields by using the BF method, the NSDP method, and the HFA method.
FrequencykBF-CPUNSDP-CPUHFA-CPU204.42203.45650.00605010.70373.24570.003810022.49324.12140.023630089.22353.76410.004440091.21434.21070.0051Figure 7
Comparisons of the NSDP method, the HFA method, and the BF method of calculating the scattered field from the patch given in Figure6 in aspects of (a) RCS values results, (b) accuracy (relative error), and (c) consumed CPU time.
(a)(b)(c)In summary, the proposed NSDP method can calculate the PO scattered fields from the quadratic concave surface in the frequency-independent computational workload and error-controllable accuracy.


## 5.1. Calculating the Assembled Concave Quadratic Triangular Patches

We consider a concave quadratic patch (see Figure4):(47)fx,y,z:z=1+0.6x2+xy+y2,with four vertex points: V1=0.8932,−13.3333, V2=22.7671,−16.6667, V3=9.8803,3.3333, and V4=−11.2201,16.6667. The incident wave direction ri=0.5,0.5,−2/2 with the amplitude Ei=−0.5,0.5,0 are considered when we treat a bistatic situation. The far field observation point lies along the direction r¯s=2/4,6/4,2/2T. The quadrilateral domain V1V2V3V4 in Figure 4(b) is decomposed into two triangular domains, ΔV1V2V3 and ΔV1V3V4, respectively.Figure 4
The original quadratic patch. (a) 3D illustration. (b) 2D illustration.
(a)(b)On invoking the affine transformation, we obtain the transformed new quadratic patch. The new vertex points for the transformed new patch areV˜1=−1.9331,−2.9957,V˜2=3.9531,−3.8366,V˜3=3.1122,1.2087,V˜4=−1.0923,4.5723.The concave quadratic patch can be calculated through two assembled triangular patches based on the NSDP method introduced above. The computational results and the corresponding consumed CPU time with increasing wavenumberk are also calculated by the high frequency asymptotic (HFA) method and the brute force (BF). Specifically, from Figure 5(a), it can be inferred that the computational results of the NSDP method agree well with those from BF method and the HFA method. Furthermore, Figure 5(b) demonstrates that the NSDP method is able to calculate the scattered fields from concave surfaces in a frequency-independent manner.Figure 5
Results and CPU time of the scattered fields with the increasing frequencyk. (a) Comparisons of the RCS values produced by the NSDP method, the BF method, and the HFA method. (b) Comparisons of the CPU time of calculating the PO scattered fields by the NSDP method, the BF method, and the HFA method.
(a)(b)


## 5.2. The Efficiency of the NSDP Method

We consider another concave PEC patch:(48)fx,y,z:z=1+0.06x2+xy+y2,to benchmark the efficiencies of our proposed NSDP method on calculating the PO scattered fields.The projection of the surface on thexy plane is a quadrilateral domain, as shown in Figure 6, with four vertex points: V1=9.0069,−13.3791, V2=22.7671,−16.6667, V3=9.8803,3.3333, and V4=−11.2201,16,6667. After the affine transformation, the quadrilateral domain in Figure 6(a) is transformed into another quadrilateral domain in Figure 6(b). We set the parameters as follows: the wave frequency k∈0,400 and the incident wave propagates along r^i=0.5,0.5,−2/2T direction. Besides the NSDP method, both the BF method and the HFA method are applied to calculate the PO scattered fields to illustrate the efficiency of the NSDP method.Figure 6
(a) The projection of the concave surface on thexy plane: V1V2V3V4. (b) The projection VV1VV2VV3VV4 after the affine transformation.
(a)(b)Through the BF method verification, we could compare the results generated by the NSDP method and the HFA method. The computational results of the PO scattered field and the relative error generated by the NSDP and the HFA methods are shown in Table1 and Table 2, respectively. Meanwhile, the CPU time consumed by these methods is given in Table 3. Figure 7(a) illustrates the RCS values of the PO scattered field Esr produced by both methods, which match well with results from the BF method. Moreover, Figure 7(b) indicates that the accuracy can be significantly improved by around two digits 102 through the NSDP method. Figure 7(c) demonstrates the CPU time consumed by both the NSDP method and the HFA method. Clearly, the computational workload for the NSDP method is relatively low and frequency-independent.Table 1
Computational results of the RCS values (dBsm unit) of PO scattered fields from the concave patch by using the BF method, the NSDP method, and the HFA method.
FrequencykBF-RCSNSDP-RCSHFA-RCS2024.748124.756724.72945026.012226.010825.986810026.102326.110726.233930027.465627.499927.437740024.917624.924524.9925Table 2
Relative error (using the results of the BF method as standard values) of calculating the PO scattered fields by using the NSDP method and the HFA method.
FrequencykNSDPHFABF203.4478×10−44.0559×10−2—505.1694×10−59.7626×10−4—1003.2366×10−45.0415×10−3—3008.8057×10−45.4585×10−3—4006.3527×10−55.0082×10−3—Table 3
Consumed CPU time (second unit) of calculating the PO scattered fields by using the BF method, the NSDP method, and the HFA method.
FrequencykBF-CPUNSDP-CPUHFA-CPU204.42203.45650.00605010.70373.24570.003810022.49324.12140.023630089.22353.76410.004440091.21434.21070.0051Figure 7
Comparisons of the NSDP method, the HFA method, and the BF method of calculating the scattered field from the patch given in Figure6 in aspects of (a) RCS values results, (b) accuracy (relative error), and (c) consumed CPU time.
(a)(b)(c)In summary, the proposed NSDP method can calculate the PO scattered fields from the quadratic concave surface in the frequency-independent computational workload and error-controllable accuracy.


## 6. Conclusion

In this work, the NSDP method is proposed to compute the highly oscillatory PO scattered fields from concave quadratic surfaces, including both the monostatic and the bistatic cases. The considered concave surface is discretized into several triangular patches. The Lagrange interpolation polynomial approximation and the affine transformation are used to approximate the amplitude function and phase function of the PO type integral. By adopting the NSDP method to obtain corresponding NSDPs, highly oscillatory PO integrand defined on the assembled triangular patches is converted to exponentially decay integrand. Numerical examples indicate the comparison of the RCS values, the CPU time, and the relative error by the NSDP method, the HFA method, and the BF method. We make the conclusion that the proposed NSDP method for calculating the PO scattered fields is error-controllable and frequency-independent for the concave quadratic patches.



---

*Source: 5296754-2020-07-23.xml*