# Creeping Ray Tracing Algorithm for Arbitrary NURBS Surfaces Based on Adaptive Variable Step Euler Method

Authors: Song Fu; Yun-Hua Zhang; Si-Yuan He; Xi Chen; Guo-Qiang Zhu
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

Although the uniform theory of diffraction (UTD) could be theoretically applied to arbitrarilyshaped convex objects modeled by nonuniform rational B-splines (NURBS), one of the great challenges in calculation of the UTD surface diffracted fields is the difficulty in determining the geodesic paths along which the creeping waves propagate on arbitrarilyshaped NURBS surfaces. In differential geometry, geodesic paths satisfy geodesic differential equation (GDE). Hence, in this paper, a general and efficient adaptive variable step Euler method is introduced for solving the GDE on arbitrarilyshaped NURBS surfaces. In contrast with conventional Euler method, the proposed method employs a shape factor (SF)ξ to efficiently enhance the accuracy of tracing and extends the application of UTD for practical engineering. The validity and usefulness of the algorithm can be verified by the numerical results.


---

## Body

## 1. Introduction

In the high-frequency electromagnetic problems, the UTD analysis of the diffraction of electromagnetic (EM) waves [1] by smooth convex surfaces is of great importance, such as the EM compatibility, the analysis of coupling between two antennas [2], and the estimation of scattering properties. And, as known to all, ray tracing of creeping waves plays an important role in the process of obtaining the solution of UTD surface diffracted flied. Therefore it is essential to investigate the creeping ray tracing (geodesic paths) on surfaces firstly. In fact, however, except some typically geometric objects which have analytical solution of the geodesic path, it is a great challenge in determining the geodesic paths along which the creeping waves propagate on arbitrarily shaped smooth convex surfaces.Jha et al. [3–5] put forward a novel technique known as the geodesic constant method, an analytical method, and it is applicable to the creeping ray tracing on general paraboloids of revolution and hybrid surfaces of revolution. Usually, for some slightly more complex models such as simple aircraft, satellite, they are usually approximated by some typically geometric objects, such as boards, cylinders, cones, and spheres, which already have analytical solution in tracing the creeping rays, for engineering application. But it is difficult to use such typically geometric objects to approximate arbitrarily shaped objects, which apparently limits the application of UTD method.Hence it is very essential to introduce the numerical ray tracing methods of creeping waves to handle arbitrarily shaped convex objects. In [6], the models are characterized by many discrete triangular meshes. And creeping ray tracing numerically on discrete triangular meshes is presented by Surazhsky et al. However, such methods cannot be used directly for the UTD solution. Except modeling arbitrarily shaped convex object by discrete triangular meshes, it can also be described in terms of free-form parametric surfaces such as NURBS surfaces [7]. Moreover, due to the advantages of high precision in modeling with low number of patches and great flexibility in simulation algorithms, NURBS has been introduced into the area of high-frequency electromagnetic analysis [8, 9]. When object is described by NURBS surfaces, geodesic path can be obtained by solving the geodesic differential equation [10]. Hence, some numerical schemes were employed to handle the GDE. In [11], to obtain the tracks of creeping rays on NURBS surfaces, the geodesic differential equations are solved by Euler method which is the first-order scheme for differential equations. And in [12], the Runge-Kutta method of order 4 was applied to solving the GDE for determining the tracing on NURBS surfaces. Euler method is highly efficient but with low precision. On the contrary, the Runge-Kutta method is of high precision but time-consuming in the process of creeping ray tracing.In this paper, aiming at improving the accuracy and efficiency of creeping ray tracing on arbitrarily shaped NURBS surface, we present a creeping ray tracing algorithm based on a novel adaptive variable step Euler method. Since the adaptive variable step Euler method is based on the conventional Euler method, the efficiency is ensured. And in the process of numerically iteratively solving the GDE, due to the introduction of shape factorξ, the discrete step sizes are adaptively corrected timely. Hence, compared with the conventional Euler method, the proposed method can easily ensure the accuracy of creeping ray tracing on arbitrarily shaped NURBS surface. Namely, it is more suitable for engineering applications. The efficiency and accuracy of the proposed method will be analyzed and verified by comparing the results of different methods.


## 2. NURBS Modeling for Arbitrarily Shaped Convex Objects

NURBS are generations of nonrational B-splines and rational and nonrational Bezier curves and surfaces [10]. And NURBS surface is the rational generalization of the nonrational B-spline surface. It is defined as follows:(1)u,v∈0,12⟶r→u,v=∑i=0m∑j=0nαijPijNipuNjqv∑i=0m∑j=0nαijNipuNjqv,where αij are the weights and Pij are the control points; Nip(t) are the normalized B-spline basis functions of degree p defined recursively as(2)Ni0t=1,ti≤t≤ti+1,0,others,Nipt=t-titi+p-tiNip-1t+ti+p+1-tti+p+1-ti+1Ni+1p-1t,where ti are the so called knots forming a knot vector T=t0,t1,…,tm. The degree, number of knots, and number of control points are related by the formula m=n+p+1.In many applications of computational geometry, NURBS surfaces are often expanded in terms of rational Bezier patches by applying the Cox-De Boor transform algorithm [13] since Bezier patches are numerically more stable than NURBS surfaces, while NURBS are more efficient for the storage and representation of a model. The surface points of a rational Bezier surface are given by(3)u,v∈0,12⟶r→u,v=∑i=0m∑j=0nαijPijBimuBjnv∑i=0m∑j=0nαijBimuBjnv,where Pij are the control points and αij are the weights. Bim(u) and Bjn(v) are Bernstein polynomials of degrees m and n, respectively, which are defined for all integers m(m≠0), n(n≠0) by(4)Bimu=m!i!m-i!ui1-um-i,Bjnv=n!j!n-j!vj1-vn-j.A NURBS or Bezier surfaceS can be defined as a vector-valued mapping from two-dimensional uv parameter domain to a set of three-dimensional xyz coordinates. And they can be generally described as(5)S:r→u,v=xu,v,yu,v,zu,vT.The mapping of uv-parameter domain to three-space can be seen in Figure 1.Figure 1
Parametric surface.In Figure1(b), when we fix v=v0 and let u vary, then r→(u,v0) depends on one parameter and is, therefore, a curve called a u-parameter curve. Similarly, if we fix u=u0, then the curve r→(u0,v) is a v-parameter curve. Note that both curves pass through P:r→(u0,v0) in R3. Tangent vectors for the u-parameter and v-parameter curves are given by differentiating the components of r→(u,v) with respect to u and v, respectively. We write(6)r→u′u,v=∂r→u,v∂u=xu′u,v,yu′u,v,zu′u,vT,r→v′u,v=∂r→u,v∂v=xv′u,v,yv′u,v,zv′u,vT.In this paper, ifr→u′(u,v) and r→v′(u,v) are both constant on parametric surfaces, we call the surface uniform grid surface, otherwise nonuniform grid surface. Here some typical NURBS models are given in Figure 2. The surfaces of cylinder are uniform grid surfaces, and others are nonuniform grid surfaces.Figure 2
NURBS models.
(a)
Cylinder
(b)
Cone
(c)
SphereThe surfaces of models above are digitized by networks ofuv-isoparametric curves. And the cylinder, cone, and sphere are, respectively, constructed by 4 Bezier patches highlighted in different colors. And in the following, for convenience, different Bezier patches will be shown in the same color.


## 3. Algorithm of Creeping Ray Tracing on Arbitrary NURBS Surfaces

According to the positions of the source and the observation points, the diffraction problems of convex surfaces can be classified as three types. First, the source and the observation points are both out of the surfaces and the observation points are far away from the surfaces. This case belongs to the scattering problem of the smooth convex surfaces. Second, the source is directly on the surfaces, and the observation points are far away from the surfaces. This can be called the radiation problem of antennas. Third, the source and the observation points are both on the surfaces. This is the mutual coupling problem among antennas on the surface. Hence, the ray tracing should also be divided into three cases, as in Figures3–5. In all of these cases, the ray trajectories on the surfaces are called creeping rays, which are restrained to propagate along geodesic paths. And this paper mainly focuses on the tracing of creeping rays since it is the most difficult part in the process of tracing.Figure 3
Ray tracing of the plane wave incidence.
(a)
Side view
(b)
Top viewFigure 4
Ray tracing of source on surface.
(a)
Side view
(b)
Top viewFigure 5
Ray tracing of source and observation point on surface.
(a)
Side view
(b)
Top view


### 3.1. The Start Points and Exit Points of the Creeping Rays

As shown in Figure6, at the starting point Q→0, the grazing incident condition satisfies(7)r→u0,v0-p→0r→u0,v0-p→0·n^u0,v0=0,where r→(u0,v0) is the unknown Q→0 on the surface r→(u,v) and n^(u0,v0) is the corresponding unit normal vector of S at Q→0.Figure 6
Creeping ray on surfaceS.At the exiting pointQ→s:r⃑(ui,vi),(8)r⃑ui,vi-p⃑sr⃑ui,vi-p⃑s·n^ui,vi=0.According to (7) and (8), we can obtain a number of the effective starting and exiting points of creeping rays on NURBS surfaces.


### 3.2. Creeping Ray Tracing by Solving GDE with Adaptive Variable Step Euler Method

As we know, the paths of creeping rays on arbitrarily shaped surfaces satisfy GDE. Then the problem of obtaining creeping rays is changed into solving GDE. Generally, GDE are solved by Euler method which is a simple and fast method. However, according to study, in most of the cases conventional Euler method does not work because of its low accuracy and stability.For the nonlinear problems, the step controlled correction procedure is essentially needed. Hence, in this paper, in order to improve the accuracy of ray tracing and make sure of the efficiency, an adaptive variable step Euler method is presented to solve the GDE.Firstly, GDE are given as follows:(9)d2uds2+Γ111duds2+2Γ121dudsdvds+Γ221dvds2=0,(10)d2vds2+Γ112duds2+2Γ122dudsdvds+Γ222dvds2=0,where Γmnk are the Christofel symbols of the second kind, s is arc length along the geodesic path, and u,v are parametric coordinates of creeping rays.Then, the first and second derivatives ofu in (9) can be approximately expressed by(11)u′s=us+h-ush,(12)u′′s=u′s+ξh-u′sξh=us+h+ξh-us+ξh-us+h+usξh2,where h is the step size between two adjacent discrete points, the determination of which is of great importance. ξ is shape factor (SF), which is introduced to adaptively control every discrete step size. And the value of ξ is subject to the shape of the object; more details about ξ will be presented in Section 3.3.In numerical calculation, a number of discrete pointsui,vi can be computed to express the creeping ray, where i=0,1,2,…. So, according to discretization, (11) and (12) can be rewritten as(13)ui′=ui+1-uihi,(14)ui′′=ui+2-u~i+1-ui+1+uiξi+1,ihi2,where u~i+1=u(s+ξi+1,ihi), which is a test point. Since(15)ui′=ui+1-uihi≈u~i+1-uiξi+1,ihi,then(16)u~i+1=ui+ξi+1,i·ui+1-ui.Next, substituting (16) into (14) we obtain(17)ui′′=ui+2-ui+1-ξi+1,i·ui+1-uiξi+1,ihi2.Similarly, the first and second derivatives ofv can be(18)vi′=vi+1-vihi,vi′′=vi+2-vi+1-ξi+1,i·vi+1-viξi+1,ihi2.Finally, substituting (13), (17), and (18) into the GDE (9) and (10), we have(19)ui+2=ui+1+ξi+1,i·Δu,vi+2=vi+1+ξi+1,i·Δv,where (20)Δu=ui+1-ui-Γ111ui+1-ui2+2Γ121ui+1-uivi+1-vi+Γ221vi+1-vi2,Δv=vi+1-vi-Γ112ui+1-ui2+2Γ122ui+1-uivi+1-vi+Γ222vi+1-vi2.According to (19), to calculate (ui+2,vi+2) on the creeping path, the previous two points (ui+1,vi+1), (ui,vi) and ξi+1,i are needed. Therefore the first two points (u0,v0), (u1,v1) and ξ1,0 should be calculated at the beginning of creeping ray tracing.Here, the starting pointQ0(u0,v0) can be obtained by (7). And according to the differential geometry [14], the second point is approximated by the sum of the two tangent vectors on tangent plane at the starting point Q0 shown in Figure 7.Figure 7
The determination of the second tracing point.Hence, the second point can be expressed by(21)u1=u0+h0·r^uu0,v0·t^r→u′u0,v0,v1=v0+h0·r^vu0,v0·t^r→v′u0,v0,where h0 is initial step size and r^u(u0,v0) and r^v(u0,v0) are the unit tangent vectors for the u-parameter and v-parameter curves at Q0(u0,v0). And tangent vectors r→u′(u0,v0) and r→v′(u0,v0) for the u-parameter and v-parameter curves can be obtained by (6). t^ is the unit vector along the incident direction.The general expression of shape factorξ is derived in Section 3.3. According to the expression, ξ1,0 can be obtained.After the calculation of(u0,v0), (u1,v1) and ξ1,0, (19) can iterate step by step with the increase of i.


### 3.3. The Derivation of Shape Factorξ

As we know, the smaller the step sizeh is, the more accurate that derivative approximate expression is. However, in the process of iteratively solving the discrete points on creeping ray path, the efficiency of algorithm will decrease with the increase of the number of discrete points. More importantly, the more the discrete points, the larger the accumulative error which may lead to wrong results.Clearly, for the approximation ofui′ or vi′, we could achieve greater efficiency if we could provide more discrete points where the variables are changing rapidly to achieve accuracy, while less discrete points should be provided in regions where variables are changing slowly. Therefore it is necessary to allocate discrete points during the process of Euler approximation.However the fact is that the parametric expressions(u(s),v(s)) of the creeping ray path are unknown. We cannot directly obtain the variable change speed degree which the policy of allocating discrete points is based on. So, the shape factor ξ is introduced to determine the discrete step size. The discrete value of shape factor ξi+1,i is defined to reflect the relative change rate of parameters at two adjacent points on the creeping path. Since the creeping rays include two parametric expressions u(s) and v(s), we should consider the change of u, v simultaneously for well determining ξi+1,i. Hence,(22)ξi+1,i=ui′ui+1′,ui′ui+1′≤vi′vi+1′,vi′vi+1′,else,where ui+1′, ui′, vi+1′, and vi′ are unknown, which need to be solved out. And the relationship between the SF ξ and step size h is(23)hi+1=ξi+1,i·hi.LetC be an arc length parametric curve (creeping ray) on surface S:r→(u,v) which pass through point P as shown in Figure 8 and denote that(24)r→s=r→us,vs.Let t^ be a unit tangent vector of C at P, and let the number of the discrete points P be i.Figure 8
CurveC on the parametric surface.Tangent vector for the creeping rayC at the point P can be obtained by differentiating equation (24) with respect to the arc length(25)t^i=dr→us,vsds=r→u,iduids+r→v,idvids,where r→u,i and r→v,i can be calculated by (6). Since |t^i|=1,(26)r→u,iduids+r→v,idvids=1.However, dui/ds and dvi/ds cannot be obtained from (26).Now, we approximatedui/ds and dvi/ds at point P on C by dui/ds at point P on u-parameter curve (v=v0) and dvi/ds at point P on v-parameter curve, respectively. For P on u-parameter curve, r→v,i=0, (26) can be(27)r→u,i′duids=1.For P on v-parameter curve, r→u,i=0, (26) can be(28)r→v,i′dvids=1.So, according to (27) and (28),(29)duids=1r→u,i′,dvids=1r→v,i′.Similarly, at the i+1 point on C,(30)dui+1ds=1r→u,i+1′,dvi+1ds=1r→v,i+1′.Substituting (29) and (30) into (22), finally the shape factor expression is(31)ξi+1,i=r→u,i+1′r→u,i′,r→u,i+1′r→u,i′≤r→v,i+1′r→v,i′,r→v,i+1′r→v,i′,else.So, based on the ξi+1,i, we can obtain discrete variable step adaptively.


## 3.1. The Start Points and Exit Points of the Creeping Rays

As shown in Figure6, at the starting point Q→0, the grazing incident condition satisfies(7)r→u0,v0-p→0r→u0,v0-p→0·n^u0,v0=0,where r→(u0,v0) is the unknown Q→0 on the surface r→(u,v) and n^(u0,v0) is the corresponding unit normal vector of S at Q→0.Figure 6
Creeping ray on surfaceS.At the exiting pointQ→s:r⃑(ui,vi),(8)r⃑ui,vi-p⃑sr⃑ui,vi-p⃑s·n^ui,vi=0.According to (7) and (8), we can obtain a number of the effective starting and exiting points of creeping rays on NURBS surfaces.


## 3.2. Creeping Ray Tracing by Solving GDE with Adaptive Variable Step Euler Method

As we know, the paths of creeping rays on arbitrarily shaped surfaces satisfy GDE. Then the problem of obtaining creeping rays is changed into solving GDE. Generally, GDE are solved by Euler method which is a simple and fast method. However, according to study, in most of the cases conventional Euler method does not work because of its low accuracy and stability.For the nonlinear problems, the step controlled correction procedure is essentially needed. Hence, in this paper, in order to improve the accuracy of ray tracing and make sure of the efficiency, an adaptive variable step Euler method is presented to solve the GDE.Firstly, GDE are given as follows:(9)d2uds2+Γ111duds2+2Γ121dudsdvds+Γ221dvds2=0,(10)d2vds2+Γ112duds2+2Γ122dudsdvds+Γ222dvds2=0,where Γmnk are the Christofel symbols of the second kind, s is arc length along the geodesic path, and u,v are parametric coordinates of creeping rays.Then, the first and second derivatives ofu in (9) can be approximately expressed by(11)u′s=us+h-ush,(12)u′′s=u′s+ξh-u′sξh=us+h+ξh-us+ξh-us+h+usξh2,where h is the step size between two adjacent discrete points, the determination of which is of great importance. ξ is shape factor (SF), which is introduced to adaptively control every discrete step size. And the value of ξ is subject to the shape of the object; more details about ξ will be presented in Section 3.3.In numerical calculation, a number of discrete pointsui,vi can be computed to express the creeping ray, where i=0,1,2,…. So, according to discretization, (11) and (12) can be rewritten as(13)ui′=ui+1-uihi,(14)ui′′=ui+2-u~i+1-ui+1+uiξi+1,ihi2,where u~i+1=u(s+ξi+1,ihi), which is a test point. Since(15)ui′=ui+1-uihi≈u~i+1-uiξi+1,ihi,then(16)u~i+1=ui+ξi+1,i·ui+1-ui.Next, substituting (16) into (14) we obtain(17)ui′′=ui+2-ui+1-ξi+1,i·ui+1-uiξi+1,ihi2.Similarly, the first and second derivatives ofv can be(18)vi′=vi+1-vihi,vi′′=vi+2-vi+1-ξi+1,i·vi+1-viξi+1,ihi2.Finally, substituting (13), (17), and (18) into the GDE (9) and (10), we have(19)ui+2=ui+1+ξi+1,i·Δu,vi+2=vi+1+ξi+1,i·Δv,where (20)Δu=ui+1-ui-Γ111ui+1-ui2+2Γ121ui+1-uivi+1-vi+Γ221vi+1-vi2,Δv=vi+1-vi-Γ112ui+1-ui2+2Γ122ui+1-uivi+1-vi+Γ222vi+1-vi2.According to (19), to calculate (ui+2,vi+2) on the creeping path, the previous two points (ui+1,vi+1), (ui,vi) and ξi+1,i are needed. Therefore the first two points (u0,v0), (u1,v1) and ξ1,0 should be calculated at the beginning of creeping ray tracing.Here, the starting pointQ0(u0,v0) can be obtained by (7). And according to the differential geometry [14], the second point is approximated by the sum of the two tangent vectors on tangent plane at the starting point Q0 shown in Figure 7.Figure 7
The determination of the second tracing point.Hence, the second point can be expressed by(21)u1=u0+h0·r^uu0,v0·t^r→u′u0,v0,v1=v0+h0·r^vu0,v0·t^r→v′u0,v0,where h0 is initial step size and r^u(u0,v0) and r^v(u0,v0) are the unit tangent vectors for the u-parameter and v-parameter curves at Q0(u0,v0). And tangent vectors r→u′(u0,v0) and r→v′(u0,v0) for the u-parameter and v-parameter curves can be obtained by (6). t^ is the unit vector along the incident direction.The general expression of shape factorξ is derived in Section 3.3. According to the expression, ξ1,0 can be obtained.After the calculation of(u0,v0), (u1,v1) and ξ1,0, (19) can iterate step by step with the increase of i.


## 3.3. The Derivation of Shape Factorξ

As we know, the smaller the step sizeh is, the more accurate that derivative approximate expression is. However, in the process of iteratively solving the discrete points on creeping ray path, the efficiency of algorithm will decrease with the increase of the number of discrete points. More importantly, the more the discrete points, the larger the accumulative error which may lead to wrong results.Clearly, for the approximation ofui′ or vi′, we could achieve greater efficiency if we could provide more discrete points where the variables are changing rapidly to achieve accuracy, while less discrete points should be provided in regions where variables are changing slowly. Therefore it is necessary to allocate discrete points during the process of Euler approximation.However the fact is that the parametric expressions(u(s),v(s)) of the creeping ray path are unknown. We cannot directly obtain the variable change speed degree which the policy of allocating discrete points is based on. So, the shape factor ξ is introduced to determine the discrete step size. The discrete value of shape factor ξi+1,i is defined to reflect the relative change rate of parameters at two adjacent points on the creeping path. Since the creeping rays include two parametric expressions u(s) and v(s), we should consider the change of u, v simultaneously for well determining ξi+1,i. Hence,(22)ξi+1,i=ui′ui+1′,ui′ui+1′≤vi′vi+1′,vi′vi+1′,else,where ui+1′, ui′, vi+1′, and vi′ are unknown, which need to be solved out. And the relationship between the SF ξ and step size h is(23)hi+1=ξi+1,i·hi.LetC be an arc length parametric curve (creeping ray) on surface S:r→(u,v) which pass through point P as shown in Figure 8 and denote that(24)r→s=r→us,vs.Let t^ be a unit tangent vector of C at P, and let the number of the discrete points P be i.Figure 8
CurveC on the parametric surface.Tangent vector for the creeping rayC at the point P can be obtained by differentiating equation (24) with respect to the arc length(25)t^i=dr→us,vsds=r→u,iduids+r→v,idvids,where r→u,i and r→v,i can be calculated by (6). Since |t^i|=1,(26)r→u,iduids+r→v,idvids=1.However, dui/ds and dvi/ds cannot be obtained from (26).Now, we approximatedui/ds and dvi/ds at point P on C by dui/ds at point P on u-parameter curve (v=v0) and dvi/ds at point P on v-parameter curve, respectively. For P on u-parameter curve, r→v,i=0, (26) can be(27)r→u,i′duids=1.For P on v-parameter curve, r→u,i=0, (26) can be(28)r→v,i′dvids=1.So, according to (27) and (28),(29)duids=1r→u,i′,dvids=1r→v,i′.Similarly, at the i+1 point on C,(30)dui+1ds=1r→u,i+1′,dvi+1ds=1r→v,i+1′.Substituting (29) and (30) into (22), finally the shape factor expression is(31)ξi+1,i=r→u,i+1′r→u,i′,r→u,i+1′r→u,i′≤r→v,i+1′r→v,i′,r→v,i+1′r→v,i′,else.So, based on the ξi+1,i, we can obtain discrete variable step adaptively.


## 4. Results and Discussions

### 4.1. Creeping Ray Tracing Algorithm Validation

The creeping rays can be calculated theoretically on some typical objects, such as cylinder, cone, and sphere. Therefore, the proposed method can be verified from the analytical results of these canonical objects.(1) Creeping Ray Tracing on Cylinder (Uniform Grid Surfaces). The cylinder which has a radius of 1 m and height of 3 m is considered for creeping ray tracing. We consider three different incident wave directions, all of which incident start from the same point Q(1,0,0). In order to validate the accuracy of the proposed creeping ray tracing method, it is compared against the analytical method and the conventional Euler method. Figure 9 shows the creeping ray tracing results obtained from the above mentioned methods, for three different incident directions, respectively. It is shown that, for each incident direction, the creeping ray obtained from the proposed method overlays with those from the other two methods. To further validate the proposed method, the creeping ray ending point and the length of the creeping rays are compared, as shown in Table 1.Table 1
Comparison of creeping rays on cylinder obtained from three different methods.
Source pointQ
(1, 0, 0)
Creeping ray number
1
2
3
Last point
Analytical method
(0.4651, −0.8853, 3.0)
(−0.1603, 0.9871, 3.0)
(0.8636, 0.5042, 3.0)
Euler method
(0.4636, −0.8867, 2.999)
(−0.1596, 0.9896, 2.9983)
(0.8638, 0.5028, 2.999)
Proposed method
(0.4636, −0.8867, 2.999)
(−0.1596, 0.9896, 2.9983)
(0.8638, 0.5028, 2.999)
Length (m)
Analytical method
5.9999
3.4636
3.0433
Euler method
6.0094
3.4682
3.0440
Proposed method
6.0094
3.4682
3.0440Figure 9
Creeping ray tracing on cylinder.From Table1, the numerical results by the proposed method are in good agreement with the analytical method, which can confirm its rationality. And we can find that the numerical results by the proposed method and Euler method, respectively, are identical, and this is because cylinder is so special that the u(s) and v(s) of creeping ray on its surfaces are first-order linear. Namely, the discrete step sizes can be chosen equally. And in this special case (uniform grid surfaces), ξi+1,i=1, the proposed method is reduced to conventional Euler method. (Comparing to the Euler method, the advantages of the proposed method can be demonstrated in the following examples.)(2) Creeping Ray Tracing on Cone (Nonuniform Grid Surfaces). The cone which has a radius of 1 m and height of 2 m is considered for creeping ray tracing. We consider three different incident wave directions, all of which incident start from the same point  Q(1,0,0). In Figure 10, the creeping ray tracing results on a cone are shown. The numbers 1, 2, and 3 in Figure 10(b) are three creeping rays of different incident directions.Figure 10
Creeping ray tracing on cone.
(a)
Side view
(b)
Top viewFirstly, from Table2, we can find the numerical tracing results by proposed method are in good agreement with the analytical method. As a comparison, the results obtained by the Euler method are in lower accuracy. Moreover, the creeping ray results obtained by the Euler method deviate from that of the analytical method when the creeping ray trajectory is closer to cone-tip. This is clearly shown in Figure 10 and Table 2, the Euler method results deviate from the analytical results for the ray number 2, and it even falls into false results when the ray runs too close to the cone-tip as for the ray number 3. This indicates that the conventional Euler method cannot be applied to nonuniform grids surfaces while the proposed method is capable of handling the nonuniform grids surfaces. Secondly, the proposed method is as efficient as conventional Euler method.Table 2
Comparison of creeping rays on cone obtained from three different methods.
Source pointQ
(1, 0, 0)
Creeping ray number
1
2
3
Last point
Analytical method
(−0.9320, −0.3622, 0.0)
(−0.0294, −0.9994, 0.0)
(1.0, −0.0389, 0.0)
Euler method
(−0.9290, −0.3594, 0.0106)
(0.0529, −0.9868, 0.0233)
(0.4461, −0.0106, 1.1080)
Proposed method
(−0.9321, −0.3546, 0.0085)
(−0.0092, −0.9900, 0.0205)
(0.9194, −0.0015, 0.1612)
Length (m)
Analytical method
3.1602
3.8669
4.3930
Euler method
3.1540
3.9004
3.1027
Proposed method
3.1610
3.8689
4.3956
Time (s)
Euler method
1.60
2.54
1.58
Proposed method
1.72
2.35
3.50(3) Creeping Ray Tracing on Sphere (Nonuniform Grid Surfaces). The sphere which has a radius of 1 m is considered for creeping ray tracing. According to differential geometry, the geodesic paths on sphere are orthodromes. Hence, the numerical results of sphere can be validated by the theoretical results.From Figure11(a), the result by proposed method is in good agreement with the theoretical results. And from Table 3 and Figures 11(b) and 11(c), although, with the decrease of discrete step size, the error by Euler method is reduced, time consumption is increased a lot. More importantly, if the discrete step size continues to decrease, from Figure 11(d), the result is totally wrong, which means Euler method is unstable.Table 3
Results of three creeping rays tracing on sphere.
Source pointQ
(1, 0, 0)
Methods
(a) Proposed method
(b) Euler method
(c) Euler method
(d) Euler method
Initial step size (m)
0.01
0.01
0.0005
0.0001
Last point
(0.9998, 0.0, 0.0)
(0.9829, 0.1817, 0.0018)
(0.9955, 0.08925, 0.0013)
(0.3462, 0.9337, −0.0086)
Time (s)
12.85
1.53
29.67
149.46Figure 11
Creeping ray tracing on sphere.
(a)
Proposed method
(b)
Euler method
(c)
Euler method
(d)
Euler method


### 4.2. Electromagnetic Calculation with UTD on the Base of Creeping Ray Tracing Algorithm

To assess the accuracy of UTD diffraction based on the proposed creeping ray tracing algorithm, the analytical results of canonical targets are given. Figures13 and 15 show the bistatic scattering electric fields of PEC sphere and cylinder in shadow region, respectively, and the UTD solutions are compared with the analytical results.Moreover, in reality, lots of targets are low detectable. When analyzing the scattering properties of those targets especially for some constructed by smooth convex curved surfaces, the contributions of creeping waves cannot be neglected. Here, taking an antiradar stealth screen component of a radar stealth satellite as example, since the antiradar stealth screen component is constructed by some smooth convex curved surfaces (see Figure16), the contributions of creeping waves will play important roles in this case. In the following example, the monostatic RCS of the antiradar stealth screen are studied by comparing the PO + UTD solution with the numerical solution to demonstrate the contribution of creeping waves.(1) The bistatic diffracted electric fields of a PEC cylinder in shadow region are calculated by UTD. The radius of the cylinder is 1 m, and the distanceρ between the observation point and the origin is 2 m. The frequency of incident plane wave is 3 GHz. The sketch map of the diffracted contribution is given in Figure 12.Figure 12
Diffracted contribution in shadow region.Figure 13
Scattered electric field of PEC cylinder.From Figure13, the UTD results of PEC cylinder in the shadow region based on the proposed creeping ray tracing algorithm are in good agreement with the analytical results, and the good agreement confirms the validity of the creeping ray tracing algorithm.(2) The bistatic diffracted electric fields of a PEC sphere in shadow region are calculated by UTD. The radius of the sphere is 1 m, and the distanceρ between the observation point and the origin is 2 m. The frequency of incident plane wave is 3 GHz. The sketch map of the diffracted contribution is given in Figure 14.Figure 14
Diffracted contribution in shadow region.Figure 15
Scattered electric field of PEC sphere.Figure 16
Creeping ray tracing on an antiradar stealth screen.As shown in Figure15, the UTD results of PEC sphere in the shadow region based on the proposed creeping ray tracing algorithm are in good agreement with the analytical results, which denotes the validity of the creeping ray tracing algorithm.(3) Monostatic RCS of the antiradar stealth screen component of a radar stealth satellite are calculated by the PO method combined with the UTD method based on the proposed creeping ray tracing algorithm, and these results are then compared with that of the numerical method to validate the accuracy of the proposed method. The side and 3D view of the antiradar stealth screen component are shown in Figure16, where R=1.12 m, r=0.3 m, and h=3.0 m.Figure16 shows the trajectories of creeping rays on the antiradar stealth screen component. In space, the conical tip of the antiradar stealth screen is directed to the earth. Hence we focus on the analysis of scattering property about the antiradar stealth screen nearby region of the conical tip. Here, the RCS will be calculated in the range of θ=90°~180°. The frequency of incident plane wave is 3 GHz.Next, we compare the monostatic RCS of the target in the range ofθ=90°~180° by three different methods, namely, the numerical method Multilevel Fast Multipole Algorithm (MLFMA), the PO method, and the PO + UTD method. In the PO + UTD scheme, the UTD solution based on the proposed algorithm is used to obtain the diffracted contribution of creeping rays. In the following, the monostatic RCS of the antiradar stealth screen are given in Figure 17.Figure 17
Monostatic RCS of the antiradar stealth screen.From Figure17, it is observed that the contribution of creeping waves is very obvious when the incident direction θ=145°~175°. In this range, thanks to the contributions of the creeping waves, the PO + UTD results are about 5 dBsm higher than the PO results. So, the contributions of creeping waves cannot be neglected. And compared with PO, the PO + UTD solution has a better agreement with the MLFMA solution.Hence, the rationality and validity of the creeping ray tracing algorithm presented in this paper are proved.


## 4.1. Creeping Ray Tracing Algorithm Validation

The creeping rays can be calculated theoretically on some typical objects, such as cylinder, cone, and sphere. Therefore, the proposed method can be verified from the analytical results of these canonical objects.(1) Creeping Ray Tracing on Cylinder (Uniform Grid Surfaces). The cylinder which has a radius of 1 m and height of 3 m is considered for creeping ray tracing. We consider three different incident wave directions, all of which incident start from the same point Q(1,0,0). In order to validate the accuracy of the proposed creeping ray tracing method, it is compared against the analytical method and the conventional Euler method. Figure 9 shows the creeping ray tracing results obtained from the above mentioned methods, for three different incident directions, respectively. It is shown that, for each incident direction, the creeping ray obtained from the proposed method overlays with those from the other two methods. To further validate the proposed method, the creeping ray ending point and the length of the creeping rays are compared, as shown in Table 1.Table 1
Comparison of creeping rays on cylinder obtained from three different methods.
Source pointQ
(1, 0, 0)
Creeping ray number
1
2
3
Last point
Analytical method
(0.4651, −0.8853, 3.0)
(−0.1603, 0.9871, 3.0)
(0.8636, 0.5042, 3.0)
Euler method
(0.4636, −0.8867, 2.999)
(−0.1596, 0.9896, 2.9983)
(0.8638, 0.5028, 2.999)
Proposed method
(0.4636, −0.8867, 2.999)
(−0.1596, 0.9896, 2.9983)
(0.8638, 0.5028, 2.999)
Length (m)
Analytical method
5.9999
3.4636
3.0433
Euler method
6.0094
3.4682
3.0440
Proposed method
6.0094
3.4682
3.0440Figure 9
Creeping ray tracing on cylinder.From Table1, the numerical results by the proposed method are in good agreement with the analytical method, which can confirm its rationality. And we can find that the numerical results by the proposed method and Euler method, respectively, are identical, and this is because cylinder is so special that the u(s) and v(s) of creeping ray on its surfaces are first-order linear. Namely, the discrete step sizes can be chosen equally. And in this special case (uniform grid surfaces), ξi+1,i=1, the proposed method is reduced to conventional Euler method. (Comparing to the Euler method, the advantages of the proposed method can be demonstrated in the following examples.)(2) Creeping Ray Tracing on Cone (Nonuniform Grid Surfaces). The cone which has a radius of 1 m and height of 2 m is considered for creeping ray tracing. We consider three different incident wave directions, all of which incident start from the same point  Q(1,0,0). In Figure 10, the creeping ray tracing results on a cone are shown. The numbers 1, 2, and 3 in Figure 10(b) are three creeping rays of different incident directions.Figure 10
Creeping ray tracing on cone.
(a)
Side view
(b)
Top viewFirstly, from Table2, we can find the numerical tracing results by proposed method are in good agreement with the analytical method. As a comparison, the results obtained by the Euler method are in lower accuracy. Moreover, the creeping ray results obtained by the Euler method deviate from that of the analytical method when the creeping ray trajectory is closer to cone-tip. This is clearly shown in Figure 10 and Table 2, the Euler method results deviate from the analytical results for the ray number 2, and it even falls into false results when the ray runs too close to the cone-tip as for the ray number 3. This indicates that the conventional Euler method cannot be applied to nonuniform grids surfaces while the proposed method is capable of handling the nonuniform grids surfaces. Secondly, the proposed method is as efficient as conventional Euler method.Table 2
Comparison of creeping rays on cone obtained from three different methods.
Source pointQ
(1, 0, 0)
Creeping ray number
1
2
3
Last point
Analytical method
(−0.9320, −0.3622, 0.0)
(−0.0294, −0.9994, 0.0)
(1.0, −0.0389, 0.0)
Euler method
(−0.9290, −0.3594, 0.0106)
(0.0529, −0.9868, 0.0233)
(0.4461, −0.0106, 1.1080)
Proposed method
(−0.9321, −0.3546, 0.0085)
(−0.0092, −0.9900, 0.0205)
(0.9194, −0.0015, 0.1612)
Length (m)
Analytical method
3.1602
3.8669
4.3930
Euler method
3.1540
3.9004
3.1027
Proposed method
3.1610
3.8689
4.3956
Time (s)
Euler method
1.60
2.54
1.58
Proposed method
1.72
2.35
3.50(3) Creeping Ray Tracing on Sphere (Nonuniform Grid Surfaces). The sphere which has a radius of 1 m is considered for creeping ray tracing. According to differential geometry, the geodesic paths on sphere are orthodromes. Hence, the numerical results of sphere can be validated by the theoretical results.From Figure11(a), the result by proposed method is in good agreement with the theoretical results. And from Table 3 and Figures 11(b) and 11(c), although, with the decrease of discrete step size, the error by Euler method is reduced, time consumption is increased a lot. More importantly, if the discrete step size continues to decrease, from Figure 11(d), the result is totally wrong, which means Euler method is unstable.Table 3
Results of three creeping rays tracing on sphere.
Source pointQ
(1, 0, 0)
Methods
(a) Proposed method
(b) Euler method
(c) Euler method
(d) Euler method
Initial step size (m)
0.01
0.01
0.0005
0.0001
Last point
(0.9998, 0.0, 0.0)
(0.9829, 0.1817, 0.0018)
(0.9955, 0.08925, 0.0013)
(0.3462, 0.9337, −0.0086)
Time (s)
12.85
1.53
29.67
149.46Figure 11
Creeping ray tracing on sphere.
(a)
Proposed method
(b)
Euler method
(c)
Euler method
(d)
Euler method


## 4.2. Electromagnetic Calculation with UTD on the Base of Creeping Ray Tracing Algorithm

To assess the accuracy of UTD diffraction based on the proposed creeping ray tracing algorithm, the analytical results of canonical targets are given. Figures13 and 15 show the bistatic scattering electric fields of PEC sphere and cylinder in shadow region, respectively, and the UTD solutions are compared with the analytical results.Moreover, in reality, lots of targets are low detectable. When analyzing the scattering properties of those targets especially for some constructed by smooth convex curved surfaces, the contributions of creeping waves cannot be neglected. Here, taking an antiradar stealth screen component of a radar stealth satellite as example, since the antiradar stealth screen component is constructed by some smooth convex curved surfaces (see Figure16), the contributions of creeping waves will play important roles in this case. In the following example, the monostatic RCS of the antiradar stealth screen are studied by comparing the PO + UTD solution with the numerical solution to demonstrate the contribution of creeping waves.(1) The bistatic diffracted electric fields of a PEC cylinder in shadow region are calculated by UTD. The radius of the cylinder is 1 m, and the distanceρ between the observation point and the origin is 2 m. The frequency of incident plane wave is 3 GHz. The sketch map of the diffracted contribution is given in Figure 12.Figure 12
Diffracted contribution in shadow region.Figure 13
Scattered electric field of PEC cylinder.From Figure13, the UTD results of PEC cylinder in the shadow region based on the proposed creeping ray tracing algorithm are in good agreement with the analytical results, and the good agreement confirms the validity of the creeping ray tracing algorithm.(2) The bistatic diffracted electric fields of a PEC sphere in shadow region are calculated by UTD. The radius of the sphere is 1 m, and the distanceρ between the observation point and the origin is 2 m. The frequency of incident plane wave is 3 GHz. The sketch map of the diffracted contribution is given in Figure 14.Figure 14
Diffracted contribution in shadow region.Figure 15
Scattered electric field of PEC sphere.Figure 16
Creeping ray tracing on an antiradar stealth screen.As shown in Figure15, the UTD results of PEC sphere in the shadow region based on the proposed creeping ray tracing algorithm are in good agreement with the analytical results, which denotes the validity of the creeping ray tracing algorithm.(3) Monostatic RCS of the antiradar stealth screen component of a radar stealth satellite are calculated by the PO method combined with the UTD method based on the proposed creeping ray tracing algorithm, and these results are then compared with that of the numerical method to validate the accuracy of the proposed method. The side and 3D view of the antiradar stealth screen component are shown in Figure16, where R=1.12 m, r=0.3 m, and h=3.0 m.Figure16 shows the trajectories of creeping rays on the antiradar stealth screen component. In space, the conical tip of the antiradar stealth screen is directed to the earth. Hence we focus on the analysis of scattering property about the antiradar stealth screen nearby region of the conical tip. Here, the RCS will be calculated in the range of θ=90°~180°. The frequency of incident plane wave is 3 GHz.Next, we compare the monostatic RCS of the target in the range ofθ=90°~180° by three different methods, namely, the numerical method Multilevel Fast Multipole Algorithm (MLFMA), the PO method, and the PO + UTD method. In the PO + UTD scheme, the UTD solution based on the proposed algorithm is used to obtain the diffracted contribution of creeping rays. In the following, the monostatic RCS of the antiradar stealth screen are given in Figure 17.Figure 17
Monostatic RCS of the antiradar stealth screen.From Figure17, it is observed that the contribution of creeping waves is very obvious when the incident direction θ=145°~175°. In this range, thanks to the contributions of the creeping waves, the PO + UTD results are about 5 dBsm higher than the PO results. So, the contributions of creeping waves cannot be neglected. And compared with PO, the PO + UTD solution has a better agreement with the MLFMA solution.Hence, the rationality and validity of the creeping ray tracing algorithm presented in this paper are proved.


## 5. Conclusion

An accurate and efficient creeping ray tracing algorithm based on an adaptive variable step Euler method is presented in this paper. The proposed algorithm is applicable to arbitrary NURBS surface. This adaptive variable step Euler method, based on the conventional Euler method, employs a shape factor (SF)ξ to solve the GDE. The SF ξ can timely reflect the relative change rate of the creeping ray parameters at two discretely adjacent points on surfaces. And according to ξ, every adaptive step size, which is very important for the accuracy of creeping ray tracing, can be obtained in turn. Finally based on the fast and accurate creeping ray tracing, the contribution of UTD diffraction can be calculated. Bearing in mind the accuracy and efficiency, the algorithm of NURBS-UTD will be very useful in practical engineering.



---

*Source: 604861-2015-08-12.xml*