- Email: [email protected]
Analysis of diffraction of a plane wave on a grating consisting of impedance bodies of revolution
Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
Analysis of diffraction of a plane wave on a grating consisting of impedance bodies of revolution A.G. Kyurkchan n, S.A. Manenkov Moscow Technical University of Communications and Informatics, Aviamotornaya Street, 8A, 111024 Moscow, Russia
a r t i c l e i n f o
abstract
Article history: Received 29 September 2010 Received in revised form 3 January 2011 Accepted 5 January 2011 Available online 12 January 2011
The problem of wave scattering by a grating consisting of coaxial impedance bodies of revolution is solved. The efficient numerical algorithm based on the modified null field method is offered. The method is applied both to scalar and vector formulations of the problem. The numerical results are obtained for various geometries of the grating elements. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Scattering problems Periodical gratings Null field method Analytical continuation of scattered field
1. Introduction The problem of diffraction of waves on periodical gratings is of great interest in optics, radiophysics, acoustics and other areas of science and engineering [1]. At the same time this problem is difficult for research. Thus the development of effective algorithms to solve the problem is rather actual. One of such approaches is the pattern equation method (PEM) [2,3]. It has been previously established that PEM allows to develop high-speed algorithms for similar problems in the case when the distance between the elements of the grating is very small. Besides using PEM one can get the analytical (asymptotical) solution of the considered problem (see below). However there is some restriction for the scatterer geometry consisting in that the scatterer (in particular the element of a grating) should belong to the class of slightly nonconvex bodies [4]. Thus to solve the problem alternative methods are essential. In this paper we develop the numerical algorithm based on modified null field method
n
Corresponding author. Tel.: + 7 499 236 2267. E-mail address: [email protected] (A.G. Kyurkchan).
0022-4073/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2011.01.002
(MNFM), which has previously been successfully applied in [5,6]. MNFM as well as the method of current integral equations is rigorous. It permits to obtain the characteristics of the wave field with any preassigned precision. The method is applicable to the diffraction problems for scatterers with analytical borders including non-convex bodies. Besides this approach permits to develop simple and universal algorithms adaptable for wide class of boundary problems for instance scattering by gratings, diffraction on a group of bodies, scattering of eigenmodes on obstacles in dielectric waveguides, etc. The present paper considers diffraction of a plane wave on the infinite periodic grating, consisting of impedance bodies of revolution located at one axis. We study both electromagnetic and acoustic statement of the problem. As it is specified above the problem is solved by MNFM. The null field method (NFM) has been offered for the first time by Waterman [7] (see also [8]). The modification of NFM named in the literature also as a method of T-matrix is successfully used to solve a wide class of diffraction problems (see for example [9]). The basis for the NFM is some relation (see below) which is satisfied everywhere inside the scatterer [7,8]. If we require that this relation is fulfilled on some closed surface inside the scatterer the
1344
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
initial boundary problem is reduced to the integral equation of the first kind relative to unknown current distributed on the surface of the body. In papers [5,6] it has been shown that for developing of the most high-speed and stable algorithms the auxiliary surface should be constructed by means of analytical deformation of the surface of the scatterer [10]. In [5,6] this version of NFM is called the modified null field method. In solving the considered problem the periodic Green’s function is used. We calculate the Green’s function in two ways. In the case of large distance (along coordinate which is perpendicular axis of a lattice) between point of source and point of observation, one can use a series obtained by the Poison formula. In the case when the distance is small it is possible to expand the Green’s function into a series of spherical harmonics. The coefficients of this series present integrals, which depend on grating period and the angle of incidence of the plane wave, but not on the geometry of its elements. This fact allows us to reduce the computation time, because these integrals can be evaluated prior to calculation of the matrix elements of algebraic system. Note that the algorithm for computing the Green’s function is analogous to the method proposed in [11], which considered the problem of diffraction on a body in circular waveguide (see also [1,12,13]). Of some interest is the asymptotical solution of the problem with the assumption that the period of the grating tends to infinity. By the help of PEM the integrooperator equation relative to the spectral function (the pattern) of the grating is derived. This equation permits us to obtain the asymptotic formulas for the pattern of the grating for large period of the structure. We apply combine method to find the pattern of the grating with large period based on PEM and MNFM. 2. The scalar problem Let us consider the grating consisting of identical coaxial impedance bodies of revolution. We assume that the grating has the period d. Introduce the Cartesian coordinate system and direct the z-axis along the axis of the grating (see Fig. 1). Denote by S0 the surface of the central element of the grating. We suppose that the structure is irradiated by the plane wave: u0 ¼ expðikrðsin y0 sin y cos j þ cos y0 cos yÞÞ
ð1Þ
where (r,y,j) are the spherical coordinates, k is the wave number in free space, y0 the incidence angle of the plane wave (on the strength of the axial symmetry of the problem geometry we take j0 =0). The diffraction field outside the grating satisfies the homogeneous Helmholtz equation
Du1 þ k2 u1 ¼ 0
ð2Þ
The scattered field also satisfies the Floquet periodic conditions: u1 ðr, j,z þ dÞ ¼ u1 ðr, j,zÞexpðikÞ
ð3Þ
where k ¼ kd cos y0 is the Floquet parameter and (r,j,z) are the cylindrical coordinates. The diffraction field
Fig. 1. Geometry of the problem.
satisfies the radiation condition at infinity: rffiffiffiffi 1 2 ip=4 X expðivs riws zÞ u1 ðr, j,zÞ , As ðjÞ e pffiffiffiffiffiffiffiffi vs r p s ¼ 1
r-1
ð4Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 where ws ¼ ðk þ 2psÞ=d, vs ¼ k ws . The sign of square root is chosen so that its imaginary part is not positive. On the surface of each element of the grating the impedance boundary condition u¼Z
@u @n
ð5Þ
! is satisfied. Here n is the outward normal and Z is the surface impedance. Let us apply MNFM. We assume that the surface of the element of the grating is analytical. It means that this surface can be described by an equation analytically depending on variables. The following relation is valid: " Z !! # @Gð r , r uÞ ! !! Jð r uÞ Gð r , r uÞZ dsu @nu S0 8 1 [ ! ! > > r 2 R3 Dn > u1 ð r Þ, > < n ¼ 1 ð6Þ ¼ 1 [ > ! > > u0 ð! r Þ, r 2 D > n :
\
n ¼ 1
where J ¼ kð@u=@nuÞ is the unknown current on the surface S0 and Dn is the domain inside the nth element of the grating. The function G in formula (6) is the periodical
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
Green’s function [14]: !! Gð r , r uÞ ¼
1 X
G0 ðRs ÞexpðiskÞ
ð7Þ
s ¼ 1
where expðikRs Þ G0 ðRs Þ ¼ , 4pkRs qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Rs ¼ r2 þ ru2 2rru cos c þ ðzzusdÞ2 ,
ð8Þ
Since the considered problem is periodic it is reduced to determination of the unknown current only on the surface of the central element of the grating. We state the following condition of null field at the auxiliary surface S0 located inside S0 [8] " Z !! # @Gð r , r uÞ ! !! ! Jð r uÞ Gð r , r uÞZ dsu ¼ u0 ð r Þ, @nu S0 ! r 2 S0
xðyÞ ¼ rðyÞexpðiyÞ, y ¼ t þ id, d 40, t 2 ½0, p
ð10Þ
If d = 0 we have xAC where C is a contour lying in the complex plane x; this contour corresponds to that of the axial section of the scatterer and is congruent with it. As d increases, C becomes smaller, and we obtain a new contour that can be chosen as the contour of the axial section of the surface S0. Such a deformation is possible as long as x(y) remains a one to one mapping. The points at which this one to one correspondence is violated satisfy the relation [10] ruðyÞ ¼i rðyÞ
ð11Þ
Here, we only take into account the roots that lie in the complex plane x inside the contour C. Formula (11) allows us to determine the critical values of the parameter d i.e., the values at which, in the course of deformation, the contour reaches the singular points of mapping (10), for different bodies. For example, for an ellipsoid of revolution [10], we have pffiffiffiffiffiffiffiffiffiffiffi ð12Þ dmax ¼ ln 2e2 =e where e is the eccentricity of the ellipse in the crosssection of the body. In the case of a multilobe body of revolution with q lobes, i.e., when rðyÞ ¼ c þ a cosðqyÞ we obtain
g ¼ a=c
ð13Þ
y ¼ rS sin yS sin j,
z ¼ rS cos yS
ð14Þ
where
yS ¼ arg xðyÞ, rS ¼ 9xðyÞ9
ð15Þ
To solve Eq. (9) numerically we expand the unknown current J and the Green’s function into the Fourier series: 1 X
Jðt, juÞ ¼
Im ðtÞexpðimjuÞ
ð16Þ
1 1 X Sm ðr, y,ru, yuÞexpðimcÞ 4p m ¼ 1
ð17Þ
m ¼ 1
Gðr, y,ru, yu, cÞ ¼
ð9Þ
Thus we get the integral equation of the first kind relative to the current on the surface of the central element of the grating. If one supposes that the surface S0 is nonresonant then it is evident that relation (9) is fulfilled in any point inside S0 because the function in the left side of formula (9) is analytical in the domain D0 [8]. As was mentioned above, we construct the surface S0 by an analytic deformation of the surface S0. Let the equation of the surface S0 in the spherical coordinate system has the form r = r(y). Then we introduce the complex variable
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 1 þ g2 ðq2 1Þ1 1 dmax ¼ ln , q gðq þ1Þ
In general, the critical value of the parameter d is determined numerically by solving Eq. (11). Thus, initially, it is necessary to determine the maximal value of the deformation parameter d. After determining this value, we assume that the parametric equations of the surface S0 are as follows: x ¼ rS sin yS cos j,
c ¼ jju
1345
where Sm ¼
Z 2p 1 1 X expðikRs iskimcÞ dc 2p s ¼ 1 0 kRs
ð18Þ
0
and t = y . Using formulas (9), (16) and (17) one can obtain the following infinite system of integral equations: Z p Km ðt,tÞIm ðtÞ dt ¼ Bm ðtÞ, m ¼ 0, 7 1, 7 2,. . ., 0
t 2 ½0, p
ð19Þ
where @Sm ðt,tÞ 2 @nu Bm ðtÞ ¼ im Jm ðkr sin y0 Þexpðikcos y0 zÞ ð20Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 Here wðtÞ ¼ rðtÞsin t r 2 ðtÞ þ ru ðtÞ and Jm ðxÞ is the Bessel function. Km ðt,tÞ ¼
wðtÞ
Sm ðt,tÞZ
3. The electromagnetic vector problem We again consider a grating consisting of identical bodies of revolution with period d and a center element bounded by the surface S0. Assume that electromagnetic plane wave !0 ! E ¼ p0 expðikrðsin y0 siny cosðjj0 Þ þcos y0 cos yÞÞ ð21Þ is incident on the grating. Outside the grating elements the scattered field satisfies the homogeneous Maxwell equations !1
!1
!1
r E ¼ ikB H , r H ¼
ik !1 E
B
ð22Þ
pffiffiffiffiffiffi where k ¼ k0 em, k0 is the wave number in free space, e and m are the dielectric pffiffiffiffiffiffiffiffi and magnetic permeabilities of the medium, B ¼ m=e is the wave impedance. The secondary field satisfies the Floquet condition !1 !1 E ðr, j,z þ dÞ ¼ E ðr, j,zÞexpðikÞ
ð23Þ
The formulas for the magnetic field are similar. The radiation condition at infinity in electromagnetic case
1346
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
takes the form rffiffiffiffi 1 ! !1 2 ip=4 X expðivs riws zÞ , A s ðjÞ E ðx,y,zÞ e pffiffiffiffiffiffiffiffi vs r p s ¼ 1
"
e12 m ¼
and the boundary condition at the surface S0 is ! ! ! ! ! n E ¼ ZU n ð n H Þ
ð25Þ
We again apply MNFM. Present the secondary field in the form Z !1 ! ! ! !! E ð r Þ ¼ iBr r J e ð r uÞGð r , r uÞ dsu S0
Z S0
! ! !! J h ð r uÞGð r , r uÞ dsu
ð34Þ
e13 m ¼
B 2 m @Sm k sin yPm þ ru sin yu @r 4p
e21 m ¼
B 1 @2 Sm k2 ðsin y cos yuSm þcos y sin yuPmþ Þ r @y@ru 4pi
ð35Þ
"
Z ! ! !1 ! i !! J h ð r uÞGð r , r uÞ dsu H ð r Þ¼ rr B S0 Z ! ! !! J e ð r uÞGð r , r uÞ dsu þkr
"
e22 m ¼
#
B 1 @2 Sm k2 ðsin y sin yu Sm þ cos y cos yu Pmþ Þ rru @y@yu 4pi
ð37Þ
S0
! ! þ kZðrGð n u J e ÞÞ dsu
ð28Þ
According to the MNFM we set the following condition at the auxiliary surface S0: Z ! ! ! ½iBðð J e , ruGÞrG þk2 J e GÞ n
e23 m ¼
B 2 m @Sm k cos yPm þ rru sin yu @y 4p
e31 m ¼
B 2 m @Sm k sin yu Pm þ r sin y @ru 4p
ð39Þ
e32 m ¼
B 2 m @Sm k cos yuPm þ rru sin y @yu 4p
ð40Þ
e33 m ¼
B m2 k2 Pmþ Sm 4pi rru sin y sin yu
h11 m ¼
ik @P þ sin y cos yumSm Dm cos y sin yusin y sin yu m 4pr sin y @y
ð27Þ
S0
! ! ! ! ! ! where J e ¼ n H , J h ¼ E n are the unknown electric and magnetic currents on the surface S0 of the central element of the grating. Using the boundary condition (25) on the surface of the central element of the grating we can rewrite the secondary electric field outside the scatterer as follows: Z !1 ! ! ! E ð r Þ¼ ½iBðð J e , ruGÞrG þk2 J e GÞ
ð29Þ
This formula can be written in the form Z ! !! ! ! !! ! ^ ! ^ ! ½Eð r , r uÞ J e ð r uÞ þ Z Hð r , r uÞð n u J e ð r uÞÞ dsu n ! !0 ¼n E
ik @P þ sin y sin yumSm þ Dm cos y cos yu þ sin y cos yu m @y 4pr sin y
h12 m ¼
1 X ! I m ðtÞexpðimjuÞ
! J e ðt, juÞ ¼
ð31Þ
m ¼ 1
The substitution of this series into formula (30) yields ^ the Fourier series for the tensor functions E^ and H: 1 X
E^ ¼
E^ m expðimcÞ,
m ¼ 1
^ ¼ H
1 X
h13 m ¼
þ k @Pm sin y þ D m cos y 4pr sin y @y
ik @P þ cos y cos yu mSm þ Dm sin y sin yu þ r sin y sin yu m @r 4pr sin y
ik @P þ cos y sin yumSm þ Dm sin y cos yu þ r sin y cos yu m @r 4pr sin y
h22 m ¼
e11 m ¼
h23 m ¼
h31 m ¼
k @Pmþ D m þr 4pr @r
ð47Þ
þ k @Pm 1 @Pmþ @Sm 1 @Sm cos y sin y þ sin yu sin y cos yu cos y @r @r 4p r @y r @y
ð48Þ
^ m expðimcÞ, H
B @2 Sm k2 ðcos y cos yuSm þ sin y sin yuPmþ Þ 4pi @r@ru
ð46Þ
ð32Þ
h32 m ¼
þ k @Pm 1 @Pmþ @Sm 1 @Sm cos yu sin y þ sin yu cos y cos y sin y þ @r @r 4p r @y r @y
^ m are [11] The components of the functions E^ m and H "
ð44Þ
ð45Þ
m ¼ 1
c ¼ jju
ð43Þ
ð30Þ
^ are the electric and magnetic tensor Here E^ and H Green’s functions, respectively. Expand the unknown current into the Fourier series
ð41Þ
ð42Þ
h21 m ¼
S0
ð38Þ
S0
! r 2 S0
#
ð36Þ
ð26Þ
! ! ! !0 þkZðrG ð n u J e ÞÞ dsu ¼ n E ,
#
r-1 ð24Þ
kr
B 1 @2 Sm k2 ðcos y sin yuSm þsin y cos yuPmþ Þ ru @r@yu 4pi
ð49Þ
# ð33Þ
h33 m ¼
ik @Pm 1 @Pm cos y sin y 4p @r r @y
ð50Þ
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
1347
functions:
where 7 ¼ 12 ððm þ1ÞSm þ 1 7 ðm1ÞSm1 Þ, Dm
Pm7 ¼ 12ðSm þ 1 7 Sm1 Þ
G0 ðRs7 Þ ¼
ð51Þ
1 1 X ð2l þ 1Þjl ðkRÞPl ðcos YÞhð2Þ ðk dsÞPl ð 71Þ l 4pi l ¼ 0
ð56Þ Define [15,16]: ! ! ruðtÞ ! 1 ! 1 2 I m ðtÞ ¼ Im i ru þ Im ðtÞ i yu þ Im ðtÞ ðtÞ i ju rðtÞ
ð52Þ
R0p 0
11 1 Km ðt,tÞIm ðtÞ dt þ 21 1 Km ðt,tÞIm ðtÞ dt þ
Rp R0p 0
12 2 Km ðt,tÞIm ðtÞ dt ¼ B1m ðtÞ,
In Eq. (56) Pl ðxÞ are the Legendre polynomials and jl ðxÞ, hð2Þ ðxÞ are the spherical Bessel and Hankel functions, l respectively. Expansion (56) is valid on condition that sd4R. Substitute formula (56) into Eq. (55). Thus we get Q X !! Gð r , r uÞ ¼ G0 ðRs ÞexpðiskÞ s ¼ Q
22 2 Km ðt,tÞIm ðtÞ dt ¼ B2m ðtÞ,
m ¼ 0, 71, 72,. . ., t 2 ½0, p
cosY ¼ ðzzuÞ=R ð57Þ
where the prime denotes the derivative with respect to the corresponding argument. Substitute formulas (31) and (32) into Eq. (30). Then we get the following infinite system of integral equations relative to the Fourier 1 2 harmonics Im , Im : (Rp
where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ r2 þ ru2 2rru cos c þðzzuÞ2 ,
þ ð53Þ
The kernels and the right hand sides of the system are ^ m and the primary expressed in terms of functions E^ m , H field, respectively. We do not present here these expressions, because they can be found numerically from the above formulas.
1 1 X ð2l þ 1Þjl ðkRÞPl ðcos YÞWl 4pi l ¼ 0
ð58Þ
where 1 X
Wl ¼
ðhð2Þ ðkdsÞexpðiksÞ þð1Þl hð2Þ ðkdsÞexpðiksÞÞ l l
s ¼ Q þ1
ð59Þ Therefore the coefficients Sm take the form Sm ¼
4. Numerical algorithm Let us divide the interval [0,p] of the parameter t by the points tj ¼ ðp=NÞðjð1=2ÞÞ where j ¼ 1,N and replace integrals in formula (53) with Riemann sums. Further we equate the obtained expressions in the collocation points which we choose at t ¼ tj . As a result we pass from the system of integral equations to the corresponding algebraic system concerning unknown values of the current at the points yu ¼ tj on the contour of axial cross-section of the surface S0. Consider the effective algorithm for calculation of the Green’s function or more exactly for calculation of the Fourier coefficients Sm. Note that the following expression can be derived by the help of the Poison formula [14]: Sm ¼
1 ip X ð2Þ Jm ðvs r o ÞHm ðvs r 4 Þexpðiws ðzzuÞÞ kd s ¼ 1
0
ð61Þ Substitution of formula (61) in Eq. (59) yields Wl ¼ Wlþ þ Wl
Wl7 ¼ i 7 l
Q X !! G0 ðRs ÞexpðiskÞ Gð r , r uÞ ¼ s ¼ Q
ðG0 ðRsþ ÞexpðiskÞ þ G0 ðR s ÞexpðiskÞÞ
s ¼ Q þ1
ð55Þ Rs7
ð60Þ
ð62Þ
where
are the Bessel and Hankel funcwhere Jm(x) and tions, r o ¼ minðr, ruÞ, r 4 ¼ maxðr, ruÞ, and ru,zu, r,z are the cylindrical coordinates of the point of source and the point of observation. Formula (54) can be used under condition 9rru9 4 1 If 9rru9 o 1 we write the presentation for the function G in the form
1 X
1 2p
Since series (59) converges very slowly we use the integral presentation for Hankel function in order to calculate (59). This presentation is Z p=2 þ i1 hð2Þ ðk dsÞ ¼ il expðik dscos aÞPl ðcos aÞsin a da l
ð54Þ
ð2Þ Hm ðxÞ
þ
! Z 2p expðikRs imcÞ dc expðiksÞ kRs 0 s ¼ Q ! Z 2p 1 X 1 i ð2l þ 1ÞWl jl ðkRÞPl ðcos YÞexpðimcÞ dc 2p 0 l¼0 Q X
¼ R 7 s . The choice of the number Q will be where explained below. We use the addition theorem for Bessel
Z p=2 þ i1 expðiðQ þ 1Þðkd cos a 7 kÞÞ 0
1expðiðkd cos a 7 kÞÞ
Pl ðcos aÞsin a da
ð63Þ The number Q in (60) is chosen so that L/((Q+ 1)kd)r 0.7 where L is the maximal cross size of the element of the grating. As follows from Eq. (63) the values Wl do not depend on the geometry of the scatterer (the element of the grating). Let us remark on the calculation of integrals in the first sum in formula (60) (so-called Vasil’ev S-functions). As is known, the real parts of these integrals have singularities and the imaginary parts are smooth functions [14]. We approximately replace these integrals with the sums Z 2p 1 1 expðikRs imcÞ 1 NX expðikRs ðcl Þimcl Þ dc 2p 0 kRs Nl¼0 kRs ðcl Þ ð64Þ
1348
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
where cl ¼ ð2plÞ=N, l ¼ 0,1,. . .,N1. Next, we apply the standard FFT procedure to calculate the sums (64). This method is suitable for determination of the S-functions when Rs is not very small. When Rs-0, the integrand in (64) changes very rapidly in the neighborhood of the point c =0. As a result, sums in Eq. (64) slowly converge as N increases. This difficulty may be obviated by means of a combined method consisting of two algorithms. When the quantity k2 rru sin y sin yu is about unity (no greater than 10–15), the real parts of S-functions are decomposed in hypergeometric functions [14]. When the quantity k2 rru sin y sin yu is large, the standard code of integration with the weight cos mc can be applied (we use standard program of the IMSL library on Fortran). When this approach is employed, it is necessary to take into account that Z 2p Z 1 expðikRs imcÞ 1 p expðikRs Þ dc ¼ cosðmcÞ dc 2p 0 kRs kRs p 0 ð65Þ This combined method of calculating of the S-functions makes it possible to determine these quantities with a high accuracy in a comparatively short CPU time, because the slow component of the algorithm (integration with the weight cos mc) is accomplished for a small number of matrix elements. Most of the matrix is calculated with the help of the FFT, which immediately yields the set of the integrals in (64) for various m. It is seen from formula (63) that the integrands have poles corresponding to the Floquet modes of the grating. In order to integrate along a contour that does not pass near the poles, we use the analyticity of the integrands and perform integration along the contour which makes an angle of 451 with the real axis in the neighborhood of the origin at the complex plane a [17]. Since the integrands are the smooth functions in a certain neighborhood of this contour, integration is easily performed with the use of standard numerical techniques. 5. Asymptotical solution of the problem for the large period of the grating
þ
!) h 1 @ cos g h @2 cos g h @2 cos g Jju Jyu þ Jru þ @j @j @yu @j @ju sin yu sin y
expðikru cos gÞwðyuÞ dyu dju
h @ cos g h @2 cos g h @2 cos g Jju Jyu þ Jru þ @y @y @yu @y @ju sin yu
rffiffiffiffi 1 !1 ! ip 2 ip=4 X expðivs riws zÞ E F 0 ðys , jÞ e pffiffiffiffiffiffiffiffi vs r kd p s ¼ 1
expðikru cos gÞwðyuÞ dyuju
Z 2p Z p 0
0
@ Jðyu, juÞ 1Z expðikru cos gÞwðyuÞ dyu dju @nu
ð68Þ ik2 F0, y ðy, jÞ ¼ 4p
Z 2p Z p ( 0
0
e @ cos g e @2 cos g e @2 cos g Jju B J þ J þ @y ru @y@yu yu @y @ju sin yu
!
ð70Þ
g0 ða, bÞ ¼ g00 ða, bÞ 1 8p2
Z 2p Z p Z 2p Z p=2 þ i1
k2 r 2 ðyÞg^0 ðau, bu; y, jÞsin y cos au
^ ðau, bu; y, jÞsin y ikruðyÞgu 0y
exp ikrðyÞðcos gcos auÞ sin au dau dbu dy dj þ
þ
1 X j¼1
1 8p2
0
0
0
0
(
expðijkd cos y0 Þ Z 2p Z p Z 2p Z p=2 þ i1 k2 r 2 ðyÞsin y cos gu k2 rðyÞruðyÞ 0
0
i0
0
@ cos gu @y expðijkd cos auÞg0 ðau, buÞexpðikrðyÞðcos gu cos gÞÞ sin au dau dbu dy dj þ expðijkd cos y0 Þ Z 2p Z p Z 2p Z p=2 þ i1 " 1 k2 r 2 ðyÞsin y cos gu þ 2 8p 0 0 0 0 sin y
@cosgu þ expðijkd cos auÞg0þ ðau, buÞ @y )
expðikrðyÞðcos gu þ cos gÞÞ sin audaudbu dy dj
ð71Þ
where g^ 0 ða, b; y, jÞ ¼
Z 2p Z p 0
v0 ðyu, juÞexp½ikrðyuÞcosg^ dyu dju
0
ð72Þ
ð67Þ
where ys ¼ arccosðws =kÞ. The functions g0(y,j) and ! F 0 ðy, jÞ are the scattering patterns of the central element of the grating. These functions are determined by the following expressions: 1 g0 ðy, jÞ ¼ 4p
!)
and where cos g ¼ sin yu sin y cosðjjuÞ þcos yu cos y e h h h ðJrue ,Jyeu ,Jj Þ, ðJ ,J ,J Þ are the components of electric and u ru yu ju ! magnetic currents in spherical coordinates. Note that J e ! and J h are related by formula (25). Consider the asymptotical solution of the problem with d-N i.e., we obtain the approximate formulas for the scattering pattern of the central element of the grating for large period. We restrict ourselves by the scalar problem of wave scattering on the grating consisting of ideal bodies of revolution (when Z = 0). Using pattern equation method one can derive the following integraloperator equation relative to the pattern g0(y,j):
k2 rðyÞruðyÞsin y
Of some interest is the asymptotics of the scattered field with r-N. In the cases of scalar and vector statement of the problem the following formulas are valid: rffiffiffiffi 1 ip 2 ip=4 X expðivs riws zÞ ð66Þ g0 ðys , jÞ u1 e pffiffiffiffiffiffiffiffi vs r kd p s ¼ 1
ð69Þ
! Z Z ( ik2 2p p B @ cos g e @2 cos g e @2 cos g Jje u F0, j ðy, jÞ ¼ Jru þ Jyu þ 4p 0 sin y @j @j @yu @j @ju sin yu 0
is the generalized pattern of the central element [18] and k @u @u v0 ðy, jÞ ¼ r 2 ðyÞsin y ruðyÞ sin y ð73Þ 4p @r @y r ¼ rðyÞ cosg^ ¼ fsin yu sinðjujÞsin a cos b þ½sin y cos yucos y sin yu cosðjujÞsin a sin b þ½cos y cos yu þ sin y sin yu cosðjujÞcos ag
ð74Þ
g07 ð
a, bÞ are determined by formulas The functions (72), where cos g^ should be replaced by cosg 7 ¼ ½sin y sin a cosðbjÞ 7cos y cos a
ð75Þ
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
1349
The values cosgu are determined by formulas (75) where the angles a, b should be replaced by the angles au, bu. Rewrite Eq. (71) as follows:
Twersky was the first to obtain equations of such kind (see, e.g., [19]). The fundamental difference between Eqs. (71) and (82) implies that the solution of the problem of diffraction by a grating with the use of Eq. (82) necessitates preliminary solution of the problem of scattering of a plane wave by a single element of the grating, while Eq. (71) allows immediate determination of the pattern of an element contained in a grating. However, Eq. (82) is suitable for the asymptotical solution of the problem for d-N. Applying the saddle-point method to estimate the integrals from (82) (assume that j0 =0) we obtain
g0 ða, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ
g0 ða, b; y0 ,0Þ g01 ða, b; y0 ,0Þ þS g01 ða, b; p,0Þg0 ðp,0; y0 ,0Þ
Then, g00 ða, b; y0 , j0 Þ Z 2p Z p i @ cos g0 ¼ k2 r 2 ðyÞsin y cos g0 k2 rðyÞruðyÞsin y 4p 0 @y 0 expðikrðyÞðcos g0 cos gÞÞ dy dj
ð76Þ
where cos g0 ¼ ½sin y sin y0 cosðjj0 Þ þ cos y cos y0
ð77Þ
7
Z 2p Z p Z 2p Z p=2 þ i1 2 2 1 þ k r ðyÞg^0 ðau, bu; y, jÞsin y cos au 2 8p 0 0 0 0
^ ðau, bu; y, jÞsin y exp ikrðyÞðcos aucos gÞ ikruðyÞgu 0y
sin au dau dbu dy dj
þ
1 X j¼1
(
Z 2pZ p=2 þ i1 1 g00 ða, b; pau, buÞ expðijkdcosy0 Þ 2pi 0 0
ð78Þ
Z 2p Z p Z 2p Z p=2 þ i1 1 1 k2 r 2 ðyÞg^ 0 ðau, bu; y, jÞsin y cos au þ 2 8p 0 0 0 0
ikruðyÞg^ 1u 0y ðau, bu; y, jÞsin y
ð79Þ
This equation can be written in the form
Accordingly, we have for Eq. (78)
j¼1
¼
1 lnð1expðikdð1 8 cos y0 ÞÞÞ kd
ð84Þ
g01 ðp,0; y0 ,0ÞS g01 ð0,0; p,0Þ þð1S g01 ðp,0; p,0ÞÞg01 ð0,0; y0 ,0Þ D
g0 ðp,0, y0 ,0Þ g 1 ð0,0; y0 ,0ÞS þ g01 ðp,0; 0,0Þ þð1S þ g01 ð0,0; 0,0ÞÞg01 ðp,0; y0 ,0Þ 0 D
ð85Þ
D ¼ ð1S g01 ðp,0; p,0ÞÞð1S þ g01 ð0,0; 0,0ÞÞ S S þ g01 ð0,0; p,0Þg01 ðp,0; 0,0Þ:
ð86Þ
Formulas (83)–(86) yield the explicit solution of the problem for kd b1 when the pattern of a single element of the grating is known. In particular the pattern can be found by MNFM (see below).
6. Numerical results
(
Z 2p Z p=2 þ i1 1 g00 ða, b; pau, buÞ expðijkd cos y0 Þ 2pi 0 0
g0 ðpau, bu; y0 , j0 Þexpðijkd cos auÞsin au dau dbu Z 2p Z p=2 þ i1 1 g00 ða, b; au, buÞ þ expðijkd cos y0 Þ 2pi 0 0 ) g0 ðau, bu; y0 , j0 Þexpðijkdcos auÞsin au dau dbu
s¼1
skd
ð80Þ
L g0 ða, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ þ 1 X
1 X expðiskd 7 iskÞ
where
g01 ða, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ
a, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ
S7 ¼
Here we indicate the dependence of the pattern on the incidence angles of the plane wave. For the single body (at d-N) we have
L½g01 ð
where
g0 ð0,0; y0 ,0Þ
)
exp½ikrðyÞðcos aucos gÞsin au dau dbudy dj
ð83Þ
From the Eq. (83) we get
g0 ðpau, bu; y0 , j0 Þexpðijkd cos auÞsin au dau dbu Z 2p Z p=2 þ i1 1 g00 ða, b; au, buÞ þ expðijkd cos y0 Þ 2pi 0 0 g0 ðau, bu; y0 , j0 ÞexpðijkdcosauÞsin au dau dbu
þ S þ g01 ða, b; 0,0Þg0 ð0,0; y0 ,0Þ
ð81Þ
In order to test the algorithm, we calculate the absolute value of the residual D of condition (30) as a function of the parameter t of the axial section contour of the auxiliary surface S0. As an example, we consider the problem of diffraction by a grating formed by perfectly conducting spheroids with the semiaxes ka =10, kc= 2.5. The period of the grating is kd =5.1, i.e., the spheroids are rather closely spaced. The electromagnetic problem is considered, and the incident field has the form
Because of the linearity of the operator L, we get g0 ða, b; y0( , j0 Þ ¼ g01 ða, b; y0 , j0 Þ Z 2pZ p=2 þ i1 1 X 1 þ g01 ða, b; pau, buÞ expðijkdcos y0 Þ 2pi 0 0 j¼1 g0 ðpau, bu; y0 , j0 Þexpðijkd cos auÞsin au dau dbu Z 2p Z p=2 þ i1 1 g01 ða, b; au, buÞ þ expðijkd cos y0 Þ 2pi 0 0 ) g0 ðau, bu; y0 , j0 Þexpðijkdcos auÞsin au dau dbu
ð82Þ
!0 ! E ¼ i z expðikxÞ,
!0 1! H ¼ i y expðikxÞ
B
ð87Þ
The number of discrete sources and the number of Fourier harmonics are N = 120 and M =15, respectively. The plot of the residual for two observation angles (j = 0 and j ¼ 903 ) is shown in Fig. 2. It follows from the figure that the maximal residual level does not exceed 2.5 10 3. To test the method we also consider the problem of wave scattering by the grating consisting of closely
1350
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
The distributions of the absolute value of the zero mode scattering pattern (i.e., the quantity p=ðkdÞ ! 9 F 0 ðy0 , jÞ9) for the gratings from elongated spheroids,
Fig. 2. The distribution of the residual at the contour of the axial section of the auxiliary surface S0. Scattering by the grating consisting from spheroids. The observation angle j =0 (solid curve) and j ¼ 903 (dashed curve).
Fig. 4. The angle dependence of the zero mode pattern for the scattering on the grating from elongated ideal spheroids.
Fig. 3. The angle dependence of the zero mode pattern of the grating from the closely spaced superellipsoids of revolution (solid curve) and the dependence of the pattern of the infinite cylinder (dashed curve).
spaced superellipsoids of revolution. The axial cross-section of the superellipsoid is defined by the equation x2l a
þ
z2l c
¼1
Fig. 5. The angle dependence of the zero mode pattern for the scattering on the grating from elongated spheroids with Z= B.
ð88Þ
At high magnitude of the parameter l and the small distance between the scatterers the problem of such geometry is little different from two-dimensional problem of the scattering by the infinite round cylinder (the normal incidence of the plane wave is assumed). As well known this problem has analytical solution. In Fig. 3 the dependence of ! the value 9 F 0 ðy0 , jÞ9 (curve 1) for the grating consisting of closely spaced superellipsoids is presented. The parameters of the problem are the following: ka=2.5, kc=5, l=10 and period of the grating kd=10.1. The grating is assumed to be irradiated by the plane wave (87). Curve 2 in the figure corresponds to the case of the scattering by the infinite round cylinder with the radius ka=2.5. One can see rather small difference between both dependences.
Fig. 6. The angle dependence of the zero mode pattern for the scattering on the grating from ideal spheres.
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
Fig. 7. The angle dependence of the zero mode pattern for the scattering on the grating from spheres with Z= B.
Fig. 8. The angle dependence of the zero mode pattern for the scattering on the grating from ideal double-lobes of revolution.
1351
dependences for the lattice consisting from double-lobes of revolution. The impedance Z of the grating elements is equal to zero (Figs. 4, 6 and 8) or equal to B (Figs. 5, 7 and 9). Thus we consider the grating which is formed from ideal elements or so-called black bodies. The parameters of the problem are as follows: the radius of the sphere ka=10, the spheroid semiaxes ka=2.5, kc=10, the parameters of the double-lobe of revolution q=2 and kc=6.25, ka=3.75. The period of the grating kd=20.1 for the case of spheres (spheroids) and kd=8.9 for the case of double-lobes of revolution. The incident plane wave is of form (87). Note that the solid curves in the figures demonstrate the dependences of the scattering pattern of the central element of the grating and the dashed curves correspond to the case of isolated body which is located in the free space and has the same sizes. As is seen from the figures the interaction between the grating elements is manifested to a greater extent when the impedance of the bodies is equal to zero. The maximal difference between the solid and dashed curves is observed in the case of the grating formed from the double-lobes of revolution and minimal difference corresponds to the case of the lattice from the spheroids (if the impedance of the spheroids is equal to B the solid and dashed curves practically coincide). We now turn to the results of asymptotical solution of the problem with kd b1. Figs. 10 and 11 demonstrate this solution for scalar problem of wave scattering by the grating consisting of ideal spheres (Fig. 10) and ideal spheroids (Fig. 11). In the figures the absolute value of scalar pattern of zero mode of the grating is plotted. The radius of the sphere is ka = 1, the semiaxes of the spheroid are kc= 1, ka= 0.25. The curves 1 in the figures correspond to the case of kd = 2.1 and curves 2 illustrate the pattern for kd = 4. The solid lines in the figures illustrate the solution of the problem obtained by MNFM and dotted lines demonstrate the asymptotical solution. As is seen from the figures, even for the shortest distance between the elements of the grating equal to 0.1, both methods of the solution of the problem differ by less than 5% for the grating of spheres and by approximately 10% for the grating of spheroids. When this distance increases to a value of two, the curves practically coincide.
Fig. 9. The angle dependence of the zero mode pattern for the scattering on the grating from double-lobes of revolution with Z = B.
spheres and double-lobes of revolution are displayed in Figs. 4–9. Figs. 4 and 5 correspond to the grating consisting of spheroids, Figs. 6 and 7 correspond to the grating formed from spheres and Figs. 8, 9 illustrate the angle
Fig. 10. The angle dependence of the zero mode pattern for the scattering on the grating from ideal spheres: comparison with the asymptotical solution (scalar problem).
1352
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
References
Fig. 11. The angle dependence of the zero mode pattern for the scattering on the grating from ideal spheroids: comparison with the asymptotical solution (scalar problem).
7. Conclusions Using MNFM the three-dimensional problem of wave scattering by the grating consisting of coaxial bodes of revolution is solved. The system of integral equations relative to the current on the central element of the grating is derived. The efficient algorithm for calculation of the Green’s function is offered. To test the method the residual at the auxiliary surface inside the central element is calculated. The obtained results demonstrate that the condition of null field is fulfilled with a high accuracy. The problem of wave scattering by closely spaced superellipsoids of revolution is also considered. It is shown that the pattern of this structure coincides very closely with the pattern of infinite cylinder of the corresponding radius. The dependences of the pattern of zero mode of the grating for various geometry of its elements are obtained. The asymptotical solution of the scalar problem for the large period of the structure is derived. This solution is illustrated by an example of the grating consisting of ideal spheres or spheroids. The results show good agreement between the asymptotical and rigorous solutions.
Acknowledgment This study was supported by the Russian Foundation for Basic Research, project nos. 09-02-00126 and 08-02-00621.
[1] Yasumoto K, editor. Electromagnetic theory and applications for photonic crystals. Taylor & Francis Group; 2006. [2] Kyurkchan AG. On solving the problem of plane wave diffraction by gratings. Phys—Dokl 1996;49:618–23. [3] Kyurkchan AG, Soloveichik AL. Scattering of waves by a periodic grating situated close to a planar interface. J Commun Technol Electron 2000;45:357–64. [4] Kyurkchan AG. A new integral equation in the diffraction theory. Sov Phys Dokl 1992;37:338–40. [5] Kyurkchan AG, Smirnova NI. Generalization of the method of extended boundary conditions. J Commun Technol Electron 2008;53:809–17. [6] Kyurkchan AG, Smirnova NI. Solution of wave diffraction problems by null-field method. Acoust Phys 2009;55:691–7. [7] Waterman PC. Matrix formulation of electromagnetic scattering. Proc IEEE 1965;53:805–12. [8] Colton DL, Kress R. In: Integral equation methods in scattering theory. New York: Wiley; 1983. [9] Mishchenko MI, Videen G, Babenko VA, Khlebtsov NG, Wriedt T. T-matrix theory of electromagnetic scattering by particles and its applications: a comprehensive reference database. JQSRT 2004;88: 357–406. [10] Kyurkchan AG, Minaev SA, Soloveitchik AL. A modification of the method of discrete sources based on a prior information about the singularities of the diffracted field. J Commun Technol Electron 2001;46:615–21. [11] Manenkov SA. Diffraction of a mode of a circular dielectric waveguide by a compact obstacle located inside the waveguide. J Commun Technol Electron 2008;53:747–57. [12] Yasumoto K, Yoshitomi K. Efficient calculation of lattice sums for free-space periodic Green’s function. IEEE Trans Antennas Propag 1999;47:1050–5. [13] Manenkov SA. Diffraction of an electromagnetic wave by a threedimensional planar lattice. J Commun Technol Electron 2010;55: 375–84. [14] Vasil’ev EN. In: Excitation of bodies of revolution. Moscow: Radio i Svyaz’; 1987. [15] Anyutin AP, Kyurkchan AG, Manenkov SA, Minaev SA. About 3D solution of diffraction problems by MMDS. JQSRT 2006;100:26–40. [16] Kyurkchan AG, Manenkov SA, Negorozhina ES. Simulation of wave scattering by a group of closely spaced bodies. J Commun Technol Electron 2006;51:1209–18. [17] Manenkov SA. Use of a spline approximation to solve the problem of diffraction by an unclosed screen located in a planar waveguide. J Commun Technol Electron 2007;52:1307–15. [18] Tkhang Do Dyk, Kyurkchan AG. Efficient method for solving the problems of wave diffraction by scatterers with broken boundaries. Acoust Phys 2003;49:43–50. [19] Twersky V. Multiple scattering of electromagnetic waves by arbitrary configurations. J Math Phys 1967;8:589–608.
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
Analysis of diffraction of a plane wave on a grating consisting of impedance bodies of revolution A.G. Kyurkchan n, S.A. Manenkov Moscow Technical University of Communications and Informatics, Aviamotornaya Street, 8A, 111024 Moscow, Russia
a r t i c l e i n f o
abstract
Article history: Received 29 September 2010 Received in revised form 3 January 2011 Accepted 5 January 2011 Available online 12 January 2011
The problem of wave scattering by a grating consisting of coaxial impedance bodies of revolution is solved. The efficient numerical algorithm based on the modified null field method is offered. The method is applied both to scalar and vector formulations of the problem. The numerical results are obtained for various geometries of the grating elements. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Scattering problems Periodical gratings Null field method Analytical continuation of scattered field
1. Introduction The problem of diffraction of waves on periodical gratings is of great interest in optics, radiophysics, acoustics and other areas of science and engineering [1]. At the same time this problem is difficult for research. Thus the development of effective algorithms to solve the problem is rather actual. One of such approaches is the pattern equation method (PEM) [2,3]. It has been previously established that PEM allows to develop high-speed algorithms for similar problems in the case when the distance between the elements of the grating is very small. Besides using PEM one can get the analytical (asymptotical) solution of the considered problem (see below). However there is some restriction for the scatterer geometry consisting in that the scatterer (in particular the element of a grating) should belong to the class of slightly nonconvex bodies [4]. Thus to solve the problem alternative methods are essential. In this paper we develop the numerical algorithm based on modified null field method
n
Corresponding author. Tel.: + 7 499 236 2267. E-mail address: [email protected] (A.G. Kyurkchan).
0022-4073/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2011.01.002
(MNFM), which has previously been successfully applied in [5,6]. MNFM as well as the method of current integral equations is rigorous. It permits to obtain the characteristics of the wave field with any preassigned precision. The method is applicable to the diffraction problems for scatterers with analytical borders including non-convex bodies. Besides this approach permits to develop simple and universal algorithms adaptable for wide class of boundary problems for instance scattering by gratings, diffraction on a group of bodies, scattering of eigenmodes on obstacles in dielectric waveguides, etc. The present paper considers diffraction of a plane wave on the infinite periodic grating, consisting of impedance bodies of revolution located at one axis. We study both electromagnetic and acoustic statement of the problem. As it is specified above the problem is solved by MNFM. The null field method (NFM) has been offered for the first time by Waterman [7] (see also [8]). The modification of NFM named in the literature also as a method of T-matrix is successfully used to solve a wide class of diffraction problems (see for example [9]). The basis for the NFM is some relation (see below) which is satisfied everywhere inside the scatterer [7,8]. If we require that this relation is fulfilled on some closed surface inside the scatterer the
1344
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
initial boundary problem is reduced to the integral equation of the first kind relative to unknown current distributed on the surface of the body. In papers [5,6] it has been shown that for developing of the most high-speed and stable algorithms the auxiliary surface should be constructed by means of analytical deformation of the surface of the scatterer [10]. In [5,6] this version of NFM is called the modified null field method. In solving the considered problem the periodic Green’s function is used. We calculate the Green’s function in two ways. In the case of large distance (along coordinate which is perpendicular axis of a lattice) between point of source and point of observation, one can use a series obtained by the Poison formula. In the case when the distance is small it is possible to expand the Green’s function into a series of spherical harmonics. The coefficients of this series present integrals, which depend on grating period and the angle of incidence of the plane wave, but not on the geometry of its elements. This fact allows us to reduce the computation time, because these integrals can be evaluated prior to calculation of the matrix elements of algebraic system. Note that the algorithm for computing the Green’s function is analogous to the method proposed in [11], which considered the problem of diffraction on a body in circular waveguide (see also [1,12,13]). Of some interest is the asymptotical solution of the problem with the assumption that the period of the grating tends to infinity. By the help of PEM the integrooperator equation relative to the spectral function (the pattern) of the grating is derived. This equation permits us to obtain the asymptotic formulas for the pattern of the grating for large period of the structure. We apply combine method to find the pattern of the grating with large period based on PEM and MNFM. 2. The scalar problem Let us consider the grating consisting of identical coaxial impedance bodies of revolution. We assume that the grating has the period d. Introduce the Cartesian coordinate system and direct the z-axis along the axis of the grating (see Fig. 1). Denote by S0 the surface of the central element of the grating. We suppose that the structure is irradiated by the plane wave: u0 ¼ expðikrðsin y0 sin y cos j þ cos y0 cos yÞÞ
ð1Þ
where (r,y,j) are the spherical coordinates, k is the wave number in free space, y0 the incidence angle of the plane wave (on the strength of the axial symmetry of the problem geometry we take j0 =0). The diffraction field outside the grating satisfies the homogeneous Helmholtz equation
Du1 þ k2 u1 ¼ 0
ð2Þ
The scattered field also satisfies the Floquet periodic conditions: u1 ðr, j,z þ dÞ ¼ u1 ðr, j,zÞexpðikÞ
ð3Þ
where k ¼ kd cos y0 is the Floquet parameter and (r,j,z) are the cylindrical coordinates. The diffraction field
Fig. 1. Geometry of the problem.
satisfies the radiation condition at infinity: rffiffiffiffi 1 2 ip=4 X expðivs riws zÞ u1 ðr, j,zÞ , As ðjÞ e pffiffiffiffiffiffiffiffi vs r p s ¼ 1
r-1
ð4Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 where ws ¼ ðk þ 2psÞ=d, vs ¼ k ws . The sign of square root is chosen so that its imaginary part is not positive. On the surface of each element of the grating the impedance boundary condition u¼Z
@u @n
ð5Þ
! is satisfied. Here n is the outward normal and Z is the surface impedance. Let us apply MNFM. We assume that the surface of the element of the grating is analytical. It means that this surface can be described by an equation analytically depending on variables. The following relation is valid: " Z !! # @Gð r , r uÞ ! !! Jð r uÞ Gð r , r uÞZ dsu @nu S0 8 1 [ ! ! > > r 2 R3 Dn > u1 ð r Þ, > < n ¼ 1 ð6Þ ¼ 1 [ > ! > > u0 ð! r Þ, r 2 D > n :
\
n ¼ 1
where J ¼ kð@u=@nuÞ is the unknown current on the surface S0 and Dn is the domain inside the nth element of the grating. The function G in formula (6) is the periodical
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
Green’s function [14]: !! Gð r , r uÞ ¼
1 X
G0 ðRs ÞexpðiskÞ
ð7Þ
s ¼ 1
where expðikRs Þ G0 ðRs Þ ¼ , 4pkRs qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Rs ¼ r2 þ ru2 2rru cos c þ ðzzusdÞ2 ,
ð8Þ
Since the considered problem is periodic it is reduced to determination of the unknown current only on the surface of the central element of the grating. We state the following condition of null field at the auxiliary surface S0 located inside S0 [8] " Z !! # @Gð r , r uÞ ! !! ! Jð r uÞ Gð r , r uÞZ dsu ¼ u0 ð r Þ, @nu S0 ! r 2 S0
xðyÞ ¼ rðyÞexpðiyÞ, y ¼ t þ id, d 40, t 2 ½0, p
ð10Þ
If d = 0 we have xAC where C is a contour lying in the complex plane x; this contour corresponds to that of the axial section of the scatterer and is congruent with it. As d increases, C becomes smaller, and we obtain a new contour that can be chosen as the contour of the axial section of the surface S0. Such a deformation is possible as long as x(y) remains a one to one mapping. The points at which this one to one correspondence is violated satisfy the relation [10] ruðyÞ ¼i rðyÞ
ð11Þ
Here, we only take into account the roots that lie in the complex plane x inside the contour C. Formula (11) allows us to determine the critical values of the parameter d i.e., the values at which, in the course of deformation, the contour reaches the singular points of mapping (10), for different bodies. For example, for an ellipsoid of revolution [10], we have pffiffiffiffiffiffiffiffiffiffiffi ð12Þ dmax ¼ ln 2e2 =e where e is the eccentricity of the ellipse in the crosssection of the body. In the case of a multilobe body of revolution with q lobes, i.e., when rðyÞ ¼ c þ a cosðqyÞ we obtain
g ¼ a=c
ð13Þ
y ¼ rS sin yS sin j,
z ¼ rS cos yS
ð14Þ
where
yS ¼ arg xðyÞ, rS ¼ 9xðyÞ9
ð15Þ
To solve Eq. (9) numerically we expand the unknown current J and the Green’s function into the Fourier series: 1 X
Jðt, juÞ ¼
Im ðtÞexpðimjuÞ
ð16Þ
1 1 X Sm ðr, y,ru, yuÞexpðimcÞ 4p m ¼ 1
ð17Þ
m ¼ 1
Gðr, y,ru, yu, cÞ ¼
ð9Þ
Thus we get the integral equation of the first kind relative to the current on the surface of the central element of the grating. If one supposes that the surface S0 is nonresonant then it is evident that relation (9) is fulfilled in any point inside S0 because the function in the left side of formula (9) is analytical in the domain D0 [8]. As was mentioned above, we construct the surface S0 by an analytic deformation of the surface S0. Let the equation of the surface S0 in the spherical coordinate system has the form r = r(y). Then we introduce the complex variable
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 1 þ g2 ðq2 1Þ1 1 dmax ¼ ln , q gðq þ1Þ
In general, the critical value of the parameter d is determined numerically by solving Eq. (11). Thus, initially, it is necessary to determine the maximal value of the deformation parameter d. After determining this value, we assume that the parametric equations of the surface S0 are as follows: x ¼ rS sin yS cos j,
c ¼ jju
1345
where Sm ¼
Z 2p 1 1 X expðikRs iskimcÞ dc 2p s ¼ 1 0 kRs
ð18Þ
0
and t = y . Using formulas (9), (16) and (17) one can obtain the following infinite system of integral equations: Z p Km ðt,tÞIm ðtÞ dt ¼ Bm ðtÞ, m ¼ 0, 7 1, 7 2,. . ., 0
t 2 ½0, p
ð19Þ
where @Sm ðt,tÞ 2 @nu Bm ðtÞ ¼ im Jm ðkr sin y0 Þexpðikcos y0 zÞ ð20Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 Here wðtÞ ¼ rðtÞsin t r 2 ðtÞ þ ru ðtÞ and Jm ðxÞ is the Bessel function. Km ðt,tÞ ¼
wðtÞ
Sm ðt,tÞZ
3. The electromagnetic vector problem We again consider a grating consisting of identical bodies of revolution with period d and a center element bounded by the surface S0. Assume that electromagnetic plane wave !0 ! E ¼ p0 expðikrðsin y0 siny cosðjj0 Þ þcos y0 cos yÞÞ ð21Þ is incident on the grating. Outside the grating elements the scattered field satisfies the homogeneous Maxwell equations !1
!1
!1
r E ¼ ikB H , r H ¼
ik !1 E
B
ð22Þ
pffiffiffiffiffiffi where k ¼ k0 em, k0 is the wave number in free space, e and m are the dielectric pffiffiffiffiffiffiffiffi and magnetic permeabilities of the medium, B ¼ m=e is the wave impedance. The secondary field satisfies the Floquet condition !1 !1 E ðr, j,z þ dÞ ¼ E ðr, j,zÞexpðikÞ
ð23Þ
The formulas for the magnetic field are similar. The radiation condition at infinity in electromagnetic case
1346
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
takes the form rffiffiffiffi 1 ! !1 2 ip=4 X expðivs riws zÞ , A s ðjÞ E ðx,y,zÞ e pffiffiffiffiffiffiffiffi vs r p s ¼ 1
"
e12 m ¼
and the boundary condition at the surface S0 is ! ! ! ! ! n E ¼ ZU n ð n H Þ
ð25Þ
We again apply MNFM. Present the secondary field in the form Z !1 ! ! ! !! E ð r Þ ¼ iBr r J e ð r uÞGð r , r uÞ dsu S0
Z S0
! ! !! J h ð r uÞGð r , r uÞ dsu
ð34Þ
e13 m ¼
B 2 m @Sm k sin yPm þ ru sin yu @r 4p
e21 m ¼
B 1 @2 Sm k2 ðsin y cos yuSm þcos y sin yuPmþ Þ r @y@ru 4pi
ð35Þ
"
Z ! ! !1 ! i !! J h ð r uÞGð r , r uÞ dsu H ð r Þ¼ rr B S0 Z ! ! !! J e ð r uÞGð r , r uÞ dsu þkr
"
e22 m ¼
#
B 1 @2 Sm k2 ðsin y sin yu Sm þ cos y cos yu Pmþ Þ rru @y@yu 4pi
ð37Þ
S0
! ! þ kZðrGð n u J e ÞÞ dsu
ð28Þ
According to the MNFM we set the following condition at the auxiliary surface S0: Z ! ! ! ½iBðð J e , ruGÞrG þk2 J e GÞ n
e23 m ¼
B 2 m @Sm k cos yPm þ rru sin yu @y 4p
e31 m ¼
B 2 m @Sm k sin yu Pm þ r sin y @ru 4p
ð39Þ
e32 m ¼
B 2 m @Sm k cos yuPm þ rru sin y @yu 4p
ð40Þ
e33 m ¼
B m2 k2 Pmþ Sm 4pi rru sin y sin yu
h11 m ¼
ik @P þ sin y cos yumSm Dm cos y sin yusin y sin yu m 4pr sin y @y
ð27Þ
S0
! ! ! ! ! ! where J e ¼ n H , J h ¼ E n are the unknown electric and magnetic currents on the surface S0 of the central element of the grating. Using the boundary condition (25) on the surface of the central element of the grating we can rewrite the secondary electric field outside the scatterer as follows: Z !1 ! ! ! E ð r Þ¼ ½iBðð J e , ruGÞrG þk2 J e GÞ
ð29Þ
This formula can be written in the form Z ! !! ! ! !! ! ^ ! ^ ! ½Eð r , r uÞ J e ð r uÞ þ Z Hð r , r uÞð n u J e ð r uÞÞ dsu n ! !0 ¼n E
ik @P þ sin y sin yumSm þ Dm cos y cos yu þ sin y cos yu m @y 4pr sin y
h12 m ¼
1 X ! I m ðtÞexpðimjuÞ
! J e ðt, juÞ ¼
ð31Þ
m ¼ 1
The substitution of this series into formula (30) yields ^ the Fourier series for the tensor functions E^ and H: 1 X
E^ ¼
E^ m expðimcÞ,
m ¼ 1
^ ¼ H
1 X
h13 m ¼
þ k @Pm sin y þ D m cos y 4pr sin y @y
ik @P þ cos y cos yu mSm þ Dm sin y sin yu þ r sin y sin yu m @r 4pr sin y
ik @P þ cos y sin yumSm þ Dm sin y cos yu þ r sin y cos yu m @r 4pr sin y
h22 m ¼
e11 m ¼
h23 m ¼
h31 m ¼
k @Pmþ D m þr 4pr @r
ð47Þ
þ k @Pm 1 @Pmþ @Sm 1 @Sm cos y sin y þ sin yu sin y cos yu cos y @r @r 4p r @y r @y
ð48Þ
^ m expðimcÞ, H
B @2 Sm k2 ðcos y cos yuSm þ sin y sin yuPmþ Þ 4pi @r@ru
ð46Þ
ð32Þ
h32 m ¼
þ k @Pm 1 @Pmþ @Sm 1 @Sm cos yu sin y þ sin yu cos y cos y sin y þ @r @r 4p r @y r @y
^ m are [11] The components of the functions E^ m and H "
ð44Þ
ð45Þ
m ¼ 1
c ¼ jju
ð43Þ
ð30Þ
^ are the electric and magnetic tensor Here E^ and H Green’s functions, respectively. Expand the unknown current into the Fourier series
ð41Þ
ð42Þ
h21 m ¼
S0
ð38Þ
S0
! r 2 S0
#
ð36Þ
ð26Þ
! ! ! !0 þkZðrG ð n u J e ÞÞ dsu ¼ n E ,
#
r-1 ð24Þ
kr
B 1 @2 Sm k2 ðcos y sin yuSm þsin y cos yuPmþ Þ ru @r@yu 4pi
ð49Þ
# ð33Þ
h33 m ¼
ik @Pm 1 @Pm cos y sin y 4p @r r @y
ð50Þ
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
1347
functions:
where 7 ¼ 12 ððm þ1ÞSm þ 1 7 ðm1ÞSm1 Þ, Dm
Pm7 ¼ 12ðSm þ 1 7 Sm1 Þ
G0 ðRs7 Þ ¼
ð51Þ
1 1 X ð2l þ 1Þjl ðkRÞPl ðcos YÞhð2Þ ðk dsÞPl ð 71Þ l 4pi l ¼ 0
ð56Þ Define [15,16]: ! ! ruðtÞ ! 1 ! 1 2 I m ðtÞ ¼ Im i ru þ Im ðtÞ i yu þ Im ðtÞ ðtÞ i ju rðtÞ
ð52Þ
R0p 0
11 1 Km ðt,tÞIm ðtÞ dt þ 21 1 Km ðt,tÞIm ðtÞ dt þ
Rp R0p 0
12 2 Km ðt,tÞIm ðtÞ dt ¼ B1m ðtÞ,
In Eq. (56) Pl ðxÞ are the Legendre polynomials and jl ðxÞ, hð2Þ ðxÞ are the spherical Bessel and Hankel functions, l respectively. Expansion (56) is valid on condition that sd4R. Substitute formula (56) into Eq. (55). Thus we get Q X !! Gð r , r uÞ ¼ G0 ðRs ÞexpðiskÞ s ¼ Q
22 2 Km ðt,tÞIm ðtÞ dt ¼ B2m ðtÞ,
m ¼ 0, 71, 72,. . ., t 2 ½0, p
cosY ¼ ðzzuÞ=R ð57Þ
where the prime denotes the derivative with respect to the corresponding argument. Substitute formulas (31) and (32) into Eq. (30). Then we get the following infinite system of integral equations relative to the Fourier 1 2 harmonics Im , Im : (Rp
where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ r2 þ ru2 2rru cos c þðzzuÞ2 ,
þ ð53Þ
The kernels and the right hand sides of the system are ^ m and the primary expressed in terms of functions E^ m , H field, respectively. We do not present here these expressions, because they can be found numerically from the above formulas.
1 1 X ð2l þ 1Þjl ðkRÞPl ðcos YÞWl 4pi l ¼ 0
ð58Þ
where 1 X
Wl ¼
ðhð2Þ ðkdsÞexpðiksÞ þð1Þl hð2Þ ðkdsÞexpðiksÞÞ l l
s ¼ Q þ1
ð59Þ Therefore the coefficients Sm take the form Sm ¼
4. Numerical algorithm Let us divide the interval [0,p] of the parameter t by the points tj ¼ ðp=NÞðjð1=2ÞÞ where j ¼ 1,N and replace integrals in formula (53) with Riemann sums. Further we equate the obtained expressions in the collocation points which we choose at t ¼ tj . As a result we pass from the system of integral equations to the corresponding algebraic system concerning unknown values of the current at the points yu ¼ tj on the contour of axial cross-section of the surface S0. Consider the effective algorithm for calculation of the Green’s function or more exactly for calculation of the Fourier coefficients Sm. Note that the following expression can be derived by the help of the Poison formula [14]: Sm ¼
1 ip X ð2Þ Jm ðvs r o ÞHm ðvs r 4 Þexpðiws ðzzuÞÞ kd s ¼ 1
0
ð61Þ Substitution of formula (61) in Eq. (59) yields Wl ¼ Wlþ þ Wl
Wl7 ¼ i 7 l
Q X !! G0 ðRs ÞexpðiskÞ Gð r , r uÞ ¼ s ¼ Q
ðG0 ðRsþ ÞexpðiskÞ þ G0 ðR s ÞexpðiskÞÞ
s ¼ Q þ1
ð55Þ Rs7
ð60Þ
ð62Þ
where
are the Bessel and Hankel funcwhere Jm(x) and tions, r o ¼ minðr, ruÞ, r 4 ¼ maxðr, ruÞ, and ru,zu, r,z are the cylindrical coordinates of the point of source and the point of observation. Formula (54) can be used under condition 9rru9 4 1 If 9rru9 o 1 we write the presentation for the function G in the form
1 X
1 2p
Since series (59) converges very slowly we use the integral presentation for Hankel function in order to calculate (59). This presentation is Z p=2 þ i1 hð2Þ ðk dsÞ ¼ il expðik dscos aÞPl ðcos aÞsin a da l
ð54Þ
ð2Þ Hm ðxÞ
þ
! Z 2p expðikRs imcÞ dc expðiksÞ kRs 0 s ¼ Q ! Z 2p 1 X 1 i ð2l þ 1ÞWl jl ðkRÞPl ðcos YÞexpðimcÞ dc 2p 0 l¼0 Q X
¼ R 7 s . The choice of the number Q will be where explained below. We use the addition theorem for Bessel
Z p=2 þ i1 expðiðQ þ 1Þðkd cos a 7 kÞÞ 0
1expðiðkd cos a 7 kÞÞ
Pl ðcos aÞsin a da
ð63Þ The number Q in (60) is chosen so that L/((Q+ 1)kd)r 0.7 where L is the maximal cross size of the element of the grating. As follows from Eq. (63) the values Wl do not depend on the geometry of the scatterer (the element of the grating). Let us remark on the calculation of integrals in the first sum in formula (60) (so-called Vasil’ev S-functions). As is known, the real parts of these integrals have singularities and the imaginary parts are smooth functions [14]. We approximately replace these integrals with the sums Z 2p 1 1 expðikRs imcÞ 1 NX expðikRs ðcl Þimcl Þ dc 2p 0 kRs Nl¼0 kRs ðcl Þ ð64Þ
1348
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
where cl ¼ ð2plÞ=N, l ¼ 0,1,. . .,N1. Next, we apply the standard FFT procedure to calculate the sums (64). This method is suitable for determination of the S-functions when Rs is not very small. When Rs-0, the integrand in (64) changes very rapidly in the neighborhood of the point c =0. As a result, sums in Eq. (64) slowly converge as N increases. This difficulty may be obviated by means of a combined method consisting of two algorithms. When the quantity k2 rru sin y sin yu is about unity (no greater than 10–15), the real parts of S-functions are decomposed in hypergeometric functions [14]. When the quantity k2 rru sin y sin yu is large, the standard code of integration with the weight cos mc can be applied (we use standard program of the IMSL library on Fortran). When this approach is employed, it is necessary to take into account that Z 2p Z 1 expðikRs imcÞ 1 p expðikRs Þ dc ¼ cosðmcÞ dc 2p 0 kRs kRs p 0 ð65Þ This combined method of calculating of the S-functions makes it possible to determine these quantities with a high accuracy in a comparatively short CPU time, because the slow component of the algorithm (integration with the weight cos mc) is accomplished for a small number of matrix elements. Most of the matrix is calculated with the help of the FFT, which immediately yields the set of the integrals in (64) for various m. It is seen from formula (63) that the integrands have poles corresponding to the Floquet modes of the grating. In order to integrate along a contour that does not pass near the poles, we use the analyticity of the integrands and perform integration along the contour which makes an angle of 451 with the real axis in the neighborhood of the origin at the complex plane a [17]. Since the integrands are the smooth functions in a certain neighborhood of this contour, integration is easily performed with the use of standard numerical techniques. 5. Asymptotical solution of the problem for the large period of the grating
þ
!) h 1 @ cos g h @2 cos g h @2 cos g Jju Jyu þ Jru þ @j @j @yu @j @ju sin yu sin y
expðikru cos gÞwðyuÞ dyu dju
h @ cos g h @2 cos g h @2 cos g Jju Jyu þ Jru þ @y @y @yu @y @ju sin yu
rffiffiffiffi 1 !1 ! ip 2 ip=4 X expðivs riws zÞ E F 0 ðys , jÞ e pffiffiffiffiffiffiffiffi vs r kd p s ¼ 1
expðikru cos gÞwðyuÞ dyuju
Z 2p Z p 0
0
@ Jðyu, juÞ 1Z expðikru cos gÞwðyuÞ dyu dju @nu
ð68Þ ik2 F0, y ðy, jÞ ¼ 4p
Z 2p Z p ( 0
0
e @ cos g e @2 cos g e @2 cos g Jju B J þ J þ @y ru @y@yu yu @y @ju sin yu
!
ð70Þ
g0 ða, bÞ ¼ g00 ða, bÞ 1 8p2
Z 2p Z p Z 2p Z p=2 þ i1
k2 r 2 ðyÞg^0 ðau, bu; y, jÞsin y cos au
^ ðau, bu; y, jÞsin y ikruðyÞgu 0y
exp ikrðyÞðcos gcos auÞ sin au dau dbu dy dj þ
þ
1 X j¼1
1 8p2
0
0
0
0
(
expðijkd cos y0 Þ Z 2p Z p Z 2p Z p=2 þ i1 k2 r 2 ðyÞsin y cos gu k2 rðyÞruðyÞ 0
0
i0
0
@ cos gu @y expðijkd cos auÞg0 ðau, buÞexpðikrðyÞðcos gu cos gÞÞ sin au dau dbu dy dj þ expðijkd cos y0 Þ Z 2p Z p Z 2p Z p=2 þ i1 " 1 k2 r 2 ðyÞsin y cos gu þ 2 8p 0 0 0 0 sin y
@cosgu þ expðijkd cos auÞg0þ ðau, buÞ @y )
expðikrðyÞðcos gu þ cos gÞÞ sin audaudbu dy dj
ð71Þ
where g^ 0 ða, b; y, jÞ ¼
Z 2p Z p 0
v0 ðyu, juÞexp½ikrðyuÞcosg^ dyu dju
0
ð72Þ
ð67Þ
where ys ¼ arccosðws =kÞ. The functions g0(y,j) and ! F 0 ðy, jÞ are the scattering patterns of the central element of the grating. These functions are determined by the following expressions: 1 g0 ðy, jÞ ¼ 4p
!)
and where cos g ¼ sin yu sin y cosðjjuÞ þcos yu cos y e h h h ðJrue ,Jyeu ,Jj Þ, ðJ ,J ,J Þ are the components of electric and u ru yu ju ! magnetic currents in spherical coordinates. Note that J e ! and J h are related by formula (25). Consider the asymptotical solution of the problem with d-N i.e., we obtain the approximate formulas for the scattering pattern of the central element of the grating for large period. We restrict ourselves by the scalar problem of wave scattering on the grating consisting of ideal bodies of revolution (when Z = 0). Using pattern equation method one can derive the following integraloperator equation relative to the pattern g0(y,j):
k2 rðyÞruðyÞsin y
Of some interest is the asymptotics of the scattered field with r-N. In the cases of scalar and vector statement of the problem the following formulas are valid: rffiffiffiffi 1 ip 2 ip=4 X expðivs riws zÞ ð66Þ g0 ðys , jÞ u1 e pffiffiffiffiffiffiffiffi vs r kd p s ¼ 1
ð69Þ
! Z Z ( ik2 2p p B @ cos g e @2 cos g e @2 cos g Jje u F0, j ðy, jÞ ¼ Jru þ Jyu þ 4p 0 sin y @j @j @yu @j @ju sin yu 0
is the generalized pattern of the central element [18] and k @u @u v0 ðy, jÞ ¼ r 2 ðyÞsin y ruðyÞ sin y ð73Þ 4p @r @y r ¼ rðyÞ cosg^ ¼ fsin yu sinðjujÞsin a cos b þ½sin y cos yucos y sin yu cosðjujÞsin a sin b þ½cos y cos yu þ sin y sin yu cosðjujÞcos ag
ð74Þ
g07 ð
a, bÞ are determined by formulas The functions (72), where cos g^ should be replaced by cosg 7 ¼ ½sin y sin a cosðbjÞ 7cos y cos a
ð75Þ
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
1349
The values cosgu are determined by formulas (75) where the angles a, b should be replaced by the angles au, bu. Rewrite Eq. (71) as follows:
Twersky was the first to obtain equations of such kind (see, e.g., [19]). The fundamental difference between Eqs. (71) and (82) implies that the solution of the problem of diffraction by a grating with the use of Eq. (82) necessitates preliminary solution of the problem of scattering of a plane wave by a single element of the grating, while Eq. (71) allows immediate determination of the pattern of an element contained in a grating. However, Eq. (82) is suitable for the asymptotical solution of the problem for d-N. Applying the saddle-point method to estimate the integrals from (82) (assume that j0 =0) we obtain
g0 ða, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ
g0 ða, b; y0 ,0Þ g01 ða, b; y0 ,0Þ þS g01 ða, b; p,0Þg0 ðp,0; y0 ,0Þ
Then, g00 ða, b; y0 , j0 Þ Z 2p Z p i @ cos g0 ¼ k2 r 2 ðyÞsin y cos g0 k2 rðyÞruðyÞsin y 4p 0 @y 0 expðikrðyÞðcos g0 cos gÞÞ dy dj
ð76Þ
where cos g0 ¼ ½sin y sin y0 cosðjj0 Þ þ cos y cos y0
ð77Þ
7
Z 2p Z p Z 2p Z p=2 þ i1 2 2 1 þ k r ðyÞg^0 ðau, bu; y, jÞsin y cos au 2 8p 0 0 0 0
^ ðau, bu; y, jÞsin y exp ikrðyÞðcos aucos gÞ ikruðyÞgu 0y
sin au dau dbu dy dj
þ
1 X j¼1
(
Z 2pZ p=2 þ i1 1 g00 ða, b; pau, buÞ expðijkdcosy0 Þ 2pi 0 0
ð78Þ
Z 2p Z p Z 2p Z p=2 þ i1 1 1 k2 r 2 ðyÞg^ 0 ðau, bu; y, jÞsin y cos au þ 2 8p 0 0 0 0
ikruðyÞg^ 1u 0y ðau, bu; y, jÞsin y
ð79Þ
This equation can be written in the form
Accordingly, we have for Eq. (78)
j¼1
¼
1 lnð1expðikdð1 8 cos y0 ÞÞÞ kd
ð84Þ
g01 ðp,0; y0 ,0ÞS g01 ð0,0; p,0Þ þð1S g01 ðp,0; p,0ÞÞg01 ð0,0; y0 ,0Þ D
g0 ðp,0, y0 ,0Þ g 1 ð0,0; y0 ,0ÞS þ g01 ðp,0; 0,0Þ þð1S þ g01 ð0,0; 0,0ÞÞg01 ðp,0; y0 ,0Þ 0 D
ð85Þ
D ¼ ð1S g01 ðp,0; p,0ÞÞð1S þ g01 ð0,0; 0,0ÞÞ S S þ g01 ð0,0; p,0Þg01 ðp,0; 0,0Þ:
ð86Þ
Formulas (83)–(86) yield the explicit solution of the problem for kd b1 when the pattern of a single element of the grating is known. In particular the pattern can be found by MNFM (see below).
6. Numerical results
(
Z 2p Z p=2 þ i1 1 g00 ða, b; pau, buÞ expðijkd cos y0 Þ 2pi 0 0
g0 ðpau, bu; y0 , j0 Þexpðijkd cos auÞsin au dau dbu Z 2p Z p=2 þ i1 1 g00 ða, b; au, buÞ þ expðijkd cos y0 Þ 2pi 0 0 ) g0 ðau, bu; y0 , j0 Þexpðijkdcos auÞsin au dau dbu
s¼1
skd
ð80Þ
L g0 ða, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ þ 1 X
1 X expðiskd 7 iskÞ
where
g01 ða, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ
a, b; y0 , j0 Þ ¼ g00 ða, b; y0 , j0 Þ
S7 ¼
Here we indicate the dependence of the pattern on the incidence angles of the plane wave. For the single body (at d-N) we have
L½g01 ð
where
g0 ð0,0; y0 ,0Þ
)
exp½ikrðyÞðcos aucos gÞsin au dau dbudy dj
ð83Þ
From the Eq. (83) we get
g0 ðpau, bu; y0 , j0 Þexpðijkd cos auÞsin au dau dbu Z 2p Z p=2 þ i1 1 g00 ða, b; au, buÞ þ expðijkd cos y0 Þ 2pi 0 0 g0 ðau, bu; y0 , j0 ÞexpðijkdcosauÞsin au dau dbu
þ S þ g01 ða, b; 0,0Þg0 ð0,0; y0 ,0Þ
ð81Þ
In order to test the algorithm, we calculate the absolute value of the residual D of condition (30) as a function of the parameter t of the axial section contour of the auxiliary surface S0. As an example, we consider the problem of diffraction by a grating formed by perfectly conducting spheroids with the semiaxes ka =10, kc= 2.5. The period of the grating is kd =5.1, i.e., the spheroids are rather closely spaced. The electromagnetic problem is considered, and the incident field has the form
Because of the linearity of the operator L, we get g0 ða, b; y0( , j0 Þ ¼ g01 ða, b; y0 , j0 Þ Z 2pZ p=2 þ i1 1 X 1 þ g01 ða, b; pau, buÞ expðijkdcos y0 Þ 2pi 0 0 j¼1 g0 ðpau, bu; y0 , j0 Þexpðijkd cos auÞsin au dau dbu Z 2p Z p=2 þ i1 1 g01 ða, b; au, buÞ þ expðijkd cos y0 Þ 2pi 0 0 ) g0 ðau, bu; y0 , j0 Þexpðijkdcos auÞsin au dau dbu
ð82Þ
!0 ! E ¼ i z expðikxÞ,
!0 1! H ¼ i y expðikxÞ
B
ð87Þ
The number of discrete sources and the number of Fourier harmonics are N = 120 and M =15, respectively. The plot of the residual for two observation angles (j = 0 and j ¼ 903 ) is shown in Fig. 2. It follows from the figure that the maximal residual level does not exceed 2.5 10 3. To test the method we also consider the problem of wave scattering by the grating consisting of closely
1350
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
The distributions of the absolute value of the zero mode scattering pattern (i.e., the quantity p=ðkdÞ ! 9 F 0 ðy0 , jÞ9) for the gratings from elongated spheroids,
Fig. 2. The distribution of the residual at the contour of the axial section of the auxiliary surface S0. Scattering by the grating consisting from spheroids. The observation angle j =0 (solid curve) and j ¼ 903 (dashed curve).
Fig. 4. The angle dependence of the zero mode pattern for the scattering on the grating from elongated ideal spheroids.
Fig. 3. The angle dependence of the zero mode pattern of the grating from the closely spaced superellipsoids of revolution (solid curve) and the dependence of the pattern of the infinite cylinder (dashed curve).
spaced superellipsoids of revolution. The axial cross-section of the superellipsoid is defined by the equation x2l a
þ
z2l c
¼1
Fig. 5. The angle dependence of the zero mode pattern for the scattering on the grating from elongated spheroids with Z= B.
ð88Þ
At high magnitude of the parameter l and the small distance between the scatterers the problem of such geometry is little different from two-dimensional problem of the scattering by the infinite round cylinder (the normal incidence of the plane wave is assumed). As well known this problem has analytical solution. In Fig. 3 the dependence of ! the value 9 F 0 ðy0 , jÞ9 (curve 1) for the grating consisting of closely spaced superellipsoids is presented. The parameters of the problem are the following: ka=2.5, kc=5, l=10 and period of the grating kd=10.1. The grating is assumed to be irradiated by the plane wave (87). Curve 2 in the figure corresponds to the case of the scattering by the infinite round cylinder with the radius ka=2.5. One can see rather small difference between both dependences.
Fig. 6. The angle dependence of the zero mode pattern for the scattering on the grating from ideal spheres.
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
Fig. 7. The angle dependence of the zero mode pattern for the scattering on the grating from spheres with Z= B.
Fig. 8. The angle dependence of the zero mode pattern for the scattering on the grating from ideal double-lobes of revolution.
1351
dependences for the lattice consisting from double-lobes of revolution. The impedance Z of the grating elements is equal to zero (Figs. 4, 6 and 8) or equal to B (Figs. 5, 7 and 9). Thus we consider the grating which is formed from ideal elements or so-called black bodies. The parameters of the problem are as follows: the radius of the sphere ka=10, the spheroid semiaxes ka=2.5, kc=10, the parameters of the double-lobe of revolution q=2 and kc=6.25, ka=3.75. The period of the grating kd=20.1 for the case of spheres (spheroids) and kd=8.9 for the case of double-lobes of revolution. The incident plane wave is of form (87). Note that the solid curves in the figures demonstrate the dependences of the scattering pattern of the central element of the grating and the dashed curves correspond to the case of isolated body which is located in the free space and has the same sizes. As is seen from the figures the interaction between the grating elements is manifested to a greater extent when the impedance of the bodies is equal to zero. The maximal difference between the solid and dashed curves is observed in the case of the grating formed from the double-lobes of revolution and minimal difference corresponds to the case of the lattice from the spheroids (if the impedance of the spheroids is equal to B the solid and dashed curves practically coincide). We now turn to the results of asymptotical solution of the problem with kd b1. Figs. 10 and 11 demonstrate this solution for scalar problem of wave scattering by the grating consisting of ideal spheres (Fig. 10) and ideal spheroids (Fig. 11). In the figures the absolute value of scalar pattern of zero mode of the grating is plotted. The radius of the sphere is ka = 1, the semiaxes of the spheroid are kc= 1, ka= 0.25. The curves 1 in the figures correspond to the case of kd = 2.1 and curves 2 illustrate the pattern for kd = 4. The solid lines in the figures illustrate the solution of the problem obtained by MNFM and dotted lines demonstrate the asymptotical solution. As is seen from the figures, even for the shortest distance between the elements of the grating equal to 0.1, both methods of the solution of the problem differ by less than 5% for the grating of spheres and by approximately 10% for the grating of spheroids. When this distance increases to a value of two, the curves practically coincide.
Fig. 9. The angle dependence of the zero mode pattern for the scattering on the grating from double-lobes of revolution with Z = B.
spheres and double-lobes of revolution are displayed in Figs. 4–9. Figs. 4 and 5 correspond to the grating consisting of spheroids, Figs. 6 and 7 correspond to the grating formed from spheres and Figs. 8, 9 illustrate the angle
Fig. 10. The angle dependence of the zero mode pattern for the scattering on the grating from ideal spheres: comparison with the asymptotical solution (scalar problem).
1352
A.G. Kyurkchan, S.A. Manenkov / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 1343–1352
References
Fig. 11. The angle dependence of the zero mode pattern for the scattering on the grating from ideal spheroids: comparison with the asymptotical solution (scalar problem).
7. Conclusions Using MNFM the three-dimensional problem of wave scattering by the grating consisting of coaxial bodes of revolution is solved. The system of integral equations relative to the current on the central element of the grating is derived. The efficient algorithm for calculation of the Green’s function is offered. To test the method the residual at the auxiliary surface inside the central element is calculated. The obtained results demonstrate that the condition of null field is fulfilled with a high accuracy. The problem of wave scattering by closely spaced superellipsoids of revolution is also considered. It is shown that the pattern of this structure coincides very closely with the pattern of infinite cylinder of the corresponding radius. The dependences of the pattern of zero mode of the grating for various geometry of its elements are obtained. The asymptotical solution of the scalar problem for the large period of the structure is derived. This solution is illustrated by an example of the grating consisting of ideal spheres or spheroids. The results show good agreement between the asymptotical and rigorous solutions.
Acknowledgment This study was supported by the Russian Foundation for Basic Research, project nos. 09-02-00126 and 08-02-00621.
[1] Yasumoto K, editor. Electromagnetic theory and applications for photonic crystals. Taylor & Francis Group; 2006. [2] Kyurkchan AG. On solving the problem of plane wave diffraction by gratings. Phys—Dokl 1996;49:618–23. [3] Kyurkchan AG, Soloveichik AL. Scattering of waves by a periodic grating situated close to a planar interface. J Commun Technol Electron 2000;45:357–64. [4] Kyurkchan AG. A new integral equation in the diffraction theory. Sov Phys Dokl 1992;37:338–40. [5] Kyurkchan AG, Smirnova NI. Generalization of the method of extended boundary conditions. J Commun Technol Electron 2008;53:809–17. [6] Kyurkchan AG, Smirnova NI. Solution of wave diffraction problems by null-field method. Acoust Phys 2009;55:691–7. [7] Waterman PC. Matrix formulation of electromagnetic scattering. Proc IEEE 1965;53:805–12. [8] Colton DL, Kress R. In: Integral equation methods in scattering theory. New York: Wiley; 1983. [9] Mishchenko MI, Videen G, Babenko VA, Khlebtsov NG, Wriedt T. T-matrix theory of electromagnetic scattering by particles and its applications: a comprehensive reference database. JQSRT 2004;88: 357–406. [10] Kyurkchan AG, Minaev SA, Soloveitchik AL. A modification of the method of discrete sources based on a prior information about the singularities of the diffracted field. J Commun Technol Electron 2001;46:615–21. [11] Manenkov SA. Diffraction of a mode of a circular dielectric waveguide by a compact obstacle located inside the waveguide. J Commun Technol Electron 2008;53:747–57. [12] Yasumoto K, Yoshitomi K. Efficient calculation of lattice sums for free-space periodic Green’s function. IEEE Trans Antennas Propag 1999;47:1050–5. [13] Manenkov SA. Diffraction of an electromagnetic wave by a threedimensional planar lattice. J Commun Technol Electron 2010;55: 375–84. [14] Vasil’ev EN. In: Excitation of bodies of revolution. Moscow: Radio i Svyaz’; 1987. [15] Anyutin AP, Kyurkchan AG, Manenkov SA, Minaev SA. About 3D solution of diffraction problems by MMDS. JQSRT 2006;100:26–40. [16] Kyurkchan AG, Manenkov SA, Negorozhina ES. Simulation of wave scattering by a group of closely spaced bodies. J Commun Technol Electron 2006;51:1209–18. [17] Manenkov SA. Use of a spline approximation to solve the problem of diffraction by an unclosed screen located in a planar waveguide. J Commun Technol Electron 2007;52:1307–15. [18] Tkhang Do Dyk, Kyurkchan AG. Efficient method for solving the problems of wave diffraction by scatterers with broken boundaries. Acoust Phys 2003;49:43–50. [19] Twersky V. Multiple scattering of electromagnetic waves by arbitrary configurations. J Math Phys 1967;8:589–608.