Applying time domain physical optics to acoustic wave backscattering problem

Applied Acoustics 71 (2010) 321–327

Contents lists available at ScienceDirect

Applied Acoustics journal homepage: www.elsevier.com/locate/apacoust

Applying time domain physical optics to acoustic wave backscattering problem Kookhyun Kim a, Kim Jin-Hyeong b, Dae-Seung Cho b,*, Woojae Seong c a

Dept. of Naval Architecture, Tongmyong University, 179, Sinseonno, Nam-gu, Busan 608-711, Republic of Korea Dept. of Naval Architecture and Ocean Engineering, Pusan National University, 30, Jangjeon-dong, Guemjeong-gu, Busan 609-735, Republic of Korea c Dept. of Naval Architecture and Ocean Engineering, Seoul National University, San 56-1, Sillim-dong, Kwanak-gu, Seoul 151-744, Republic of Korea b

a r t i c l e

i n f o

Article history: Received 7 March 2009 Received in revised form 30 August 2009 Accepted 26 October 2009 Available online 22 November 2009 Keywords: Transient analysis Acoustic wave backscattering Time domain physical optics Modulated Gaussian pulse Central finite difference method Hidden surface removal

a b s t r a c t The interest in transient analysis of acoustic waves has been growing in recent years, due to the advance of wide-band sonars. In this paper, a transient analysis method for acoustic backscattering signals is proposed based on the time domain physical optics (TDPO). TDPO is formulated via a theoretical inverse Fourier transform of the conventional physical optics formula used in the frequency domain wave scattering analyses. A hidden surface removal algorithm using an adaptive triangular beam method and a virtual surface concept are adopted to explain shadow effects and multiple reflections among elements, respectively. Numerical analyses are carried out for two kinds of underwater targets: a submarine pressure hull and an idealized submarine, in order to validate the proposed method. The result of the submarine pressure hull shows good agreements between the proposed method and conventional physical optics based on inverse fast Fourier transform. Additionally, the result of the idealized submarine shows that the proposed method is efficient for finding highlights including their contribution to the whole backscattering signal. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction An accurate and realistic model for transient analysis of acoustic wave backscattering is certainly necessary in developing accurate and robust signal processing algorithm, particularly for active sonar systems extensively used for detection of underwater targets such as submarines, underwater mines and torpedoes, and for geological or geodesic survey of the ocean bottom. For the purpose, a Fourier transform-based approach has been classically adopted in the transient backscattering analysis of acoustic waves, where numerical simulations in the frequency domain are performed at equally spaced frequencies, and then backscattered signals in the time domain are recovered via inverse Fourier transform [1]. As the Fourier transform-based approach is a computationally expensive task for scattering problems of acoustically large and complex targets, numerical techniques directly operated in the time domain for a more efficient analysis have recently received considerable attention. Particularly, in electromagnetics, various kinds of time domain numerical techniques such as the time domain electric field integral equation (TD-EFIE) and time domain magnetic field integral equation (TD-MFIE) [2], and finite difference time domain methods (FDTD) [3] have been applied to the problems of determining the transient response from scatterers. Main advantages of these kinds of time domain techniques are that the involved phenomena can be * Corresponding author. Tel.: +82 51 510 2482; fax: +82 51 512 8836. E-mail address: [email protected] (D.-S. Cho). 0003-682X/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.apacoust.2009.10.006

understood with a deeper physical insight and wideband analyses are performed by executing the time domain algorithm just once. Nevertheless, they have inherent disadvantage for electromagnetically large and complex targets in terms of computing cost. In order to overcome this, the time domain physical optics (TDPO) has been developed in electromagnetics [4]. TDPO has been derived by theoretically taking the inverse Fourier transform of the conventional frequency-domain physical optics (FDPO) formula used in electromagnetic wave scattering problem. Since the introduction by Sun and Rusch [4] in electromagnetics, TDPO has been used in various applications such as beam forming of radar reflector [5] and electromagnetic wave scattering analysis [6]. Recently, hybrid approaches such as FDTD–TDPO [7] and FEBI (finite element–boundary integral)–TDPO [8] have been proposed for analysis of electromagnetic scattering. The basic assumption in physical optics is that only the illuminated surfaces directly contribute toward the scattering energy. Therefore, an efficient algorithm to find the illuminated surfaces should be adopted for more accurate analysis, particularly, of complex-shaped underwater targets. For the purpose, noteworthy algorithms such as multi-resolution shooting and bouncing rays [9] and an object precision method [10] based on the adaptive triangular beam method [11] have been applied in conventional FDPO. Nevertheless, it is hard to find an example applying these algorithms in TDPO. In this paper, application of TDPO to the transient analysis of acoustic wave backscattering signals is attempted. An acoustic

322

K. Kim et al. / Applied Acoustics 71 (2010) 321–327

TDPO is newly formulated by taking the theoretical inverse Fourier transform to FDPO derived by applying the Kirchhoff approximation to the Kirchhoff–Helmholtz integral equation. Moreover, to find the surfaces illuminated by initially incident or multi-reflected acoustic waves, a hidden surface algorithm [10] and a virtual surface concept [12] are adopted. To validate the proposed method, numerical analyses are carried out for a submarine pressure hull and an idealized submarine. The results of the submarine pressure hull are compared with the results from the classically so-called inverse fast Fourier transform physical optics (IFFT-PO) method [1], where the inverse fast Fourier transform is taken with respect to the frequency responses calculated by the frequency-domain physical optics. Then, for the idealized submarine, the locations of major highlights and their contributions to the whole backscattering signal are investigated. 2. Formulation of time domain physical optics 2.1. Review of frequency-domain physical optics Consider an acoustic plane wave incident towards the origin of the coordinate system located within a target (see Fig. 1). The acoustic field scattered by the target is obtained by the Kirchhoff–Helmholtz integral equation, in the frequency domain, as [13]

Ps ð~ r; xÞ ¼

 Z  @Pð~ r0 Þ @Gð~ r;~ r0 Þ Gð~ r;~ r0 Þ  Pð~ r0 Þ dS @n @n S

ð1Þ

where Ps ð~ r; xÞ is the frequency response of the acoustic field scattered from the target at any receiver position ~ r. Pð~ r 0 Þ is the frequency response of the acoustic field induced at any position vector ~ r;~ r 0 Þ and @Gð~ r;~ r 0 Þ=@n are the Green’s funcr 0 on the target surface. Gð~ tion and its partial derivative with respect to the normal direction of the target surface, respectively. By assuming far-field condition r;~ r 0 Þ=@n are approximated as follows: (r ! 1), Gð~ r; ~ r0 Þ and @Gð~ ~

0

ejks ð~r~r Þ ejkr j~ks ~r0  e Gðr; r Þ ¼ 4pj~ r ~ r0 j 4pr ~ ~0

  @ 1 r;~ r 0 Þ  jk cos /s Gð~ r;~ r0 Þ Gð~ r;~ r 0 Þ ¼  jk þ cos /s Gð~ @n r

ð2Þ

ð3Þ

pffiffiffiffiffiffiffi where j ¼ 1 and r ¼ j~ rj is the distance between the field point ks j is the wave number of scattering and the centre of target. ks ¼ j~ ^ and scattering wave. /s is the angle between unit normal vector n direction vector ~ r ~ r 00 . By adopting the Kirchhoff approximation, the acoustic field Pð~ rÞ on the target surface illuminated by the incident wave can be expressed as

Pð~ r 0 Þ ¼ P i ð~ r 0 Þ þ Ps ð~ r0 Þ

ð4Þ

where P i ð~ r 0 Þ and Ps ð~ r 0 Þ are the incident acoustic pressure and the scattered acoustic pressure at any position on the target surface, respectively. These are approximated by assuming far-field condition as follows: ~

0

Pi ð~ r 0 Þ ¼ P0 ð~ r 0 ; xÞejki ~r @ ~ 0 r 0 Þ  jk cos /i Pi ð~ r 0 Þ  jk cos /i P0 ð~ r 0 ; xÞejki ~r Pi ð~ @n ~ 0 Ps ð~ r 0 Þ ¼ CPi ð~ rÞ ¼ CP0 ð~ r 0 ; xÞejki ~r @ @ ~ 0 r 0 Þ  C Pi ð~ r 0 Þ  jk cos /i CP 0 ð~ r0 ; xÞejki ~r Ps ð~ @n @n

ð5Þ ð6Þ ð7Þ ð8Þ

where C is the complex reflection coefficient with respect to the an^ and incident direction vecgle /i between the unit normal vector n r0  ~ r0 is the source position vector and ki ð¼ j~ ki jÞ is the wave tor ~ r0  ~ r0 ; xÞ is the frequency response of the number of incident wave. P0 ð~ acoustic incident wave at the source position and can be obtained by taking the Fourier transform of the time history of incident acoustic wave. By assuming that both acoustic source and receiver are located at the same position on the z-axis far from the target, i.e. far-field condition and mono-static case, the following relations are satiski ¼ ~ ks , ~ ks  ð~ r ~ r 0 Þ ¼ r  z0 , and ~ r0 ¼ ~ r, where z0 fied: /i ¼ /s ¼ /, ~ is the z-component of ~ r 0 . By using these relations and substituting Eqs. (2)–(4) into the right side of Eq. (1), the scattered acoustic pressure at the position ~ r is rewritten by

Ps ð~ r; xÞ ¼ 

jk jkr e 2pr

Z

0

Cejkz Pi ð~ r 0 Þ cos /dS

ð9Þ

S

where kð¼ jki j ¼ jks jÞ is the scalar wavenumber. When the target is divided by a finite number of flat plate elements as shown in Fig. 2, Eq. (9) can be approximated by the following equation:

Ps ð~ r; xÞ ¼ 

Z M X 0 jk r; xÞejkr Cm e2jkz cos /m dS P0 ð~ 2pr Sm m¼1

ð10Þ

where M is the total number of flat plate elements comprising the target. Cm and Sm are the complex reflection coefficient and the surface area of each element, respectively. /m is the angle between the unit normal vector of mth triangular element and incident direction

Fig. 1. The coordinate system employed for the incidence and scattering of acoustic waves.

Fig. 2. Target divided into a finite number of flat triangular elements.

323

K. Kim et al. / Applied Acoustics 71 (2010) 321–327

Fig. 3. Depiction of hidden surface removal algorithm. Fig. 4. Virtual surface concept.

vector ~ r ~ r 0 . The integral term is the phase integral and can be calculated analytically [12]. Fourier synthetic approach has been classically adopted in the transient backscattering analysis of acoustic waves [1], in which numerical simulations in frequency domain are performed by using Eq. (10) for an interested set of equally spaced frequencies and then backscattered signals in the time domain are recovered via the inverse fast Fourier transform (IFFT). 2.2. Time domain physical optics formulation and its numerical implementation TDPO can be easily derived by theoretically taking the inverse Fourier transform of the FDPO. By taking a the theoretical inverse the Fourier transform with respect to time to the Eq. (9), an acoustic field in the time domain is derived as follows:

Z 1 1 Ps ð~ r xÞejwt dx 2p 1  Z 1 Z 0 1 jk jkr ¼  Cejkz P i ð~ rÞ cos /dS ejwt dx e 2p 1 2pr S  Z 1 Z  0 1 1 jwðtrz c Þ ¼ C Pi ð~ rÞjwe dx cos /dS 2prc s 2p 1  Z 1 Z  1 1 @ rz0 C P  ið~ rÞ ejwðt c Þ dx cos /dS ¼ 2prc S 2p 1 @t   Z 1 Z 1 @ 1 rz0 C Pi ð~ rÞejwðt c Þ dx cos /dS ¼ 2prc S @t 2p 1    Z 1 @ r  z0 pi ~ r; t  cos /dS C ¼ 2prc S @t c

Ps ð~ r; tÞ ¼

ð11Þ

where P s ð~ r; tÞ is the time-history of the scattered acoustic wave r 0 ; tÞ is the time-history of incident acoustic wave from the target, pi ð~ at a certain position on the target surface, x is the circular frequency, t is the time, and c is the speed of acoustic wave. Looking at Eq. (11), the scattered field using TDPO simply becomes a function of the time history of acoustic incident wave arriving on the target surface. TDPO does not require time consumable efforts such as Fourier transform of acoustic incident wave, frequency responses calculation of the target and inverse Fourier transform to get time histories of acoustic scattering wave, and complex value integration, as well. In order to calculate the time derivative of acoustic incident wave, a central finite difference scheme [14] is applied as defined:

@ p ð~ r 0 ; t  DtÞ  pi ð~ r0 ; t þ DtÞ fp ð~ r 0 ; tÞg ¼ i @t i 2 Dt where Dt is the infinitesimal of time.

ð12Þ

Fig. 5. Realistic underwater targets: (a) submarine pressure hull and (b) idealized submarine.

By substituting Eq. (12) into Eq. (11), a transient acoustic wave scattering for the target divided by a finite number of flat triangular elements can be rewritten as

ps ð~ r; tÞ ¼ 

Z M X 1 Cm F m ð~ r 0 ; tÞ cos /m dS 4prcDt m¼1 Sm

ð13Þ

where

    r  z0 r  z0  pi ~ r 0 ; t  Dt  r 0 ; t þ Dt  F m ð~ r 0 ; tÞ ¼ pi ~ c c

ð14Þ

For calculation of the surface integral of each element, the coordinate transform with z0 ¼ z01 a þ z02 b þ z03 c [15] is applied, where z01 , z02 and z03 are z-values of vertexes constructing the flat triangle

324

K. Kim et al. / Applied Acoustics 71 (2010) 321–327

Fig. 6. Acoustic signal backscattered from the submarine pressure hull when rotation angle of the target are (a) 0°, (b) 30°, (c) 60°, and (d) 90°.

element. a, b and c are the axis of the transformed coordinates satisfying the relation a þ b þ c ¼ 1. Then, Eq. (13) can be rewritten as

ps ð~ r; tÞ ¼ 

Z Z M X 1 Cm F m f~ r0 ða; b; 1  a  bÞ; t gSm 4prcDt m¼1 b a

 cos /m dadb

ð15Þ

where Sm is the area of the mth triangular element. The Eq. (15) is integrated by using Gaussian quadrature [14] for each triangle element for the target. Meanwhile, it is assumed in physical optics that only the illuminated surfaces directly contribute to the whole scattering energy. Therefore an efficient algorithm to find the illuminated surfaces should be adopted for more accurate analysis, particularly for complex underwater targets. In this paper, a hidden surface algorithm [10] with an adaptive triangular beam method [11] and a virtual surface concept [12,16] are adopted to explain shadow effects and multiple reflections among elements, respectively. Fig. 3 briefly delineates the hidden surface removal algorithm. As an example for the hidden surface algorithm, consider that a triangular element A lies in front of a rectangular hexahedron and the acoustic plane wave is incident on the target along the direction indicated in the figure, where the surface of the rectangular hexahedron is composed of 12 triangular elements and A is the shadow area projected from the triangular element A on the front-side of the hexahedron. All the triangular elements except the two elements on the front-side are removed. Then, the remaining elements are split into sub-elements referring to the left and right boundary. Finally, hidden sub-elements are removed. As a result, illuminated elements only remain. In addition, Fig. 4 shows the virtual surface concept, which produces a virtual surface corresponding to a half of the total path length traced by multiple reflections by the following equation:

z0n ¼

Ln þ Dn 2

ð16Þ

where n is the vertex number; n ¼ 1; 2; 3 for triangular element. z0n is the z-component of the nth vertex of the virtual surface . Ln is the path length from 1st reflection to 2nd reflection. Dn is the path length difference in z-direction between 1st reflection and 2nd reflection of nth vertex. With virtual elements generated from Eq. (16), one can calculate time domain backscattering signals by using Eq. (11). 3. Numerical investigations In order to validate the proposed method, transient analyses are carried out for realistic underwater targets such as a pressure hull of a submarine and an idealized submarine as shown in Fig. 5. In the analysis, it is assumed that the surface of each target is rigid, i.e. C ¼ 1. A transient signal of acoustic incident waves is generated by using a modulated Gaussian pulse defined as

pi ð~ r; tÞ ¼ sinðxm tÞefðtt0 Þ=T s g

2

ð17Þ

where xm ¼ 2pfm is the modulation circular frequency, fm is the modulation frequency. t 0 and T s are the time delay and the pulse width, respectively. In this paper, as acoustic plane wave incidence is assumed, the incident wave arriving at a certain position on the target surface r0 ; tÞ can be written as: from the source pi ð~

  r  z0 r; t  pi ð~ r 0 ; tÞ ¼ pi ~ c    2 rz0 r  z0 efðtt0  c Þ=T s g ¼ sin xm t  c

ð18Þ

K. Kim et al. / Applied Acoustics 71 (2010) 321–327

325

Fig. 7. Acoustic signal backscattered from the idealized submarine when rotation angle of the target are (a) 0°, (b) 30°, (c) 60°, (d) 90°, (e) 120°, (f) 150°, and (g) 180°.

The submarine pressure hull is a model devised to verify the target strength analysis software developed in various countries [13] and is modeled with 4600 triangle elements in this paper.

The length of the model is about 44.1 m and both the width and height are 7.0 m. The modulation frequency is set to be 500 Hz, the pulse width and time delay are set to be 3 and 6 ms,

326

K. Kim et al. / Applied Acoustics 71 (2010) 321–327

respectively. The distance between the sonar and target is set to 1500 m. Fig. 6 compares the analysis results by TDPO and those by IFFT-PO [1], where rotational angles of the target are 0°, 30°, 60° and 90°. The analysis results by TDPO are coincident with those of IFFT-PO. These prove that the numerical implementation of the proposed TDPO is proper. Additionally, the following features are observed: just one pulsed signal appears for the rotation angle of 0°, since the fore-end cap contributes as a main backscatter while the aft-end cap and the cylinder body do not contribute due to the shadow effect. For the rotation angle of 30°, the fore-end cap is still a major backscattered, but the backscattered signal from the aftend cap begins to appear a little. The effect by cylinder body does not appear because the incident angle is still off-normal to that. For the rotation angle of 60°, the backscattered signal takes on a similar aspect with that for the rotation angle of 30°. While the amplitude of the signal from the fore-end cap decreases, the amplitude of the signal from the aft-end increases due to the change of the relative aspect angle. A pulse signal with large amplitude appears for the rotation angle of 90°, which is from the cylinder body, and conceals the signals from the fore-end cap and the aft-end cap. The idealized submarine is a model to simulate target strength originated by Defense Research and Development Canada (DRDC) [17] and is modeled with 9330 triangle elements in this paper. The length of the model is 62 m, and the width and height are 7.0 m and 10.5 m, respectively. The other analysis conditions are same with those of submarine pressure hull model except that three rotation angles of 120°, 150°, and 180° are added. Fig. 7 shows the analysis result of the idealized submarine for various rotation angles: 0°, 30°, 60°, 90°, 120°, 150°, and 180°, in which the x–z plane view of the objective target is transparently superimposed onto the analysis result in order to identify intuitionally the location of highlights. For the rotation angle of 0°, two signals are observed, which are backscattered from the bow structure and the sail of the target, respectively. For the rotation angle of 30°, four signals are observed. Two front signals come from the bow structure and the sail and a successive signal with relatively small amplitude is due to the double reflection between the sail and the main body so that the signal is delayed and elongated. The last signal is due to the triple reflection between the tail structure and the tail: tail structure–tail–tail structure. For the rotation angle of 60°, while the signal shows similar patterns with that for the rotation angle of 30°, the signal seems to be more compressed. This is due to the reduced characteristic length. For the rotation angle of 90°, a strong backscattered signal appears coming from the main body and the sail. For the rotation angle of 120°, two signals with large amplitude and a signal with small amplitude appear. Two front signals are from the tail structure and the sail, respectively, and the successive signal is from bow structure. For the rotation angle of 150°, two signals with large amplitude and three signals with small amplitude appear. Two front signals are from the tail structure and the sail, respectively, similarly for those of the rotation angle of 120°. Third successive signals are due to double reflections between the sail and the main body, and the rest of the signals are from the main body and the bow structure. Finally, for the rotation angle of 180°, a signal with large amplitude and two signals with small amplitude appear. The front signal is due to the tail structure, the tail and double reflection between each other. The successive signals are from the sail. As a whole, transient signals obtained from this give useful information on major highlights of the target. From the analysis results above, it is recognized that intuitive identification of the highlights including their locations and contributions can be carried out by using the proposed method. Finally, one has to investigate how much computation costs could be reduced by using the proposed method, comparing with IFFT-PO. As mentioned above, the proposed method need not com-

Table 1 Computation cost comparison between TDPO and IFFT-PO for the idealized submarine model.

a

Incident angle (°)

TDPO (s)

IFFT-POa (s)

[TDPO]/[IFFT-PO] (%)

0 30 60 90 120 150 180 Summation

23 19 11 8 12 14 43 130

79 61 38 31 34 44 142 429

29 31 29 26 35 32 30 30

Computational costs for FFT processes not included.

putational efforts on FFT for transforming time domain to frequency-domain results and IFFT for opposite direction, and on complex value integrations as well. The conventional IFFT-PO requires two (2) more processes differently with the proposed method; one (1) Fourier transformation of incident wave field and one (1) inverse Fourier transformation of frequency-domain results. In this study, an algorithm so called radix-2 FFT is used, of which the number of computation per each process is estimated to ðN=2Þlog2 N (N: number of frequency-domain data). As a result, although neglecting efforts for complex value integration of IFFTPO and required memories, one can reduce computational costs by Nlog2 N, using TDPO. Table 1 compares computational costs between both approaches for the idealized submarine model, which shows that TDPO reduce computational cost to 30% in average comparing with IFFT-PO. For information, calculations to measure computational costs have been carried out on a desktop computer (CPU = 2.0 GHz), and the computational cost of IFFT-PO is not include that of FFT processes mentioned above. 4. Concluding remarks A transient analysis method of acoustic wave backscattered from underwater targets was established based on the time domain physical optics applied to electromagnetic wave scattering problems. Moreover, a hidden surface removal algorithm and a virtual surface concept are adopted to explain the shadow effects and multiple reflections among elements. Numerical analyses showed that the proposed method is accurate and applicable to acoustic wave backscattering problems of underwater targets such as a submarine pressure hull and an idealized submarine. Particularly, it is recognized that intuitive identification of the highlights including their locations and contributions can be efficiently carried out by using the proposed method. Additionally, the proposed method is more efficient than a conventional IFFT based approach, in terms of computational cost. Acknowledgement This work was supported by Advanced Ship Engineering Research Center of the Korea Science and Engineering Foundation. References [1] Kim K, Cho DS, Seong W. Simulation of time-domain acoustic wave signals backscattered from underwater targets. J Acoust Soc Korea 2008;27:140–8 [in Korean]. [2] Jung BH, Sarkar TK, Salazar-Palma M. Time domain EFIE and MFIE formulations for analysis of transient electromagnetic scattering from 3-D dielectric objects. Progr Electromagn Res 2004;49:113–42. [3] Kunz KS, Lee KM. A three-dimensional finite-difference solution of the external response of an aircraft to a complex transient EM environment; part I – the method and its implementation and part II – comparison of predictions and measurements. IEEE Trans Electromagn Compatibl 1978;EMC-20:328–41.

K. Kim et al. / Applied Acoustics 71 (2010) 321–327 [4] Sun EY, Rusch WVT. Time-domain physical optics. IEEE Trans Anten Propag 1994;42(1):9–15. [5] Gego CG, Hasselmann JV, Moreira JS. Time-domain analysis of a reflector antenna illuminated by a Gaussian pulse. J Microw Optoelectron 1999;1(4):20–8. [6] Yang LX, Ge DB, Wei B. Time-domain physical-optics method for transient scattering analysis. In: the 2006 4th Asia-Pacific conference on environmental electromagnetics; 2006. p. 861–4. [7] Yang LX, Ge DB, Wei B. FDTD/TDPO hybrid approach for analysis of the EM scattering of combinative objects. Progr Electromagn Res 2007;76:275–84. [8] Faghihi F, Heydari H. A combination of time domain finite element-boundary integral with time domain physical optics for calculation of electromagnetic scattering of 3-D structures. Progr Electromagn Res 2008;79:463–74. [9] Suk S, Seo TI, Park HS, Kim HT. Multiresolution grid algorithm in the SBR and its application to the RCS calculation. Microw Optical Technol Lett 2001;29(6):394–7. [10] Kim K, Cho DS, Kim CJ. High frequency acoustic scattering analysis of underwater target. J Soc Naval Architects Korea 2005;42(5):528–33 [in Korean].

327

[11] Cho DS, Sung SK, Kim JH, Choi JH, Park IK. A study on the indoor sound-field analysis by adaptive triangular beam method. Trans Korean Soc Noise Vib Eng 2003;13(3):217–24 [in Korean]. [12] Klement D, Preissner J, Stein V. Special problems in applying the physical optics method for backscatter computation of complicated objects. IEEE Trans Anten Propag 1988;36(2):228–37. [13] Schneider HG, Berg G, Gilroy L. Acoustic scattering by a submarine: results from a benchmark target strength simulation workshop, ICSV10; 2003. p. 2475–82. [14] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C: the art of scientific computing. 2nd ed. Cambridge University Press; 1992. [15] Bathe KJ. Finite element procedures. Prentice-Hall: Englewood Cliffs; 1995. [16] Kim K, Cho JH, Kim JH, Cho DS. A fast estimation of sonar cross section of acoustically large and complex underwater targets using a deterministic scattering center model. Appl Acoust 2009;70(5):653–60. [17] Nell CW, Gilroy LE. An improved BASIS model for the BeTSSi submarine. DRDC Atlantic TR 2003-199; 2003.