Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Shadows created by an array of acoustically hard cylinders above a rigid plane Clark, David 1994

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1994-0549.pdf [ 2.39MB ]
JSON: 831-1.0065320.json
JSON-LD: 831-1.0065320-ld.json
RDF/XML (Pretty): 831-1.0065320-rdf.xml
RDF/JSON: 831-1.0065320-rdf.json
Turtle: 831-1.0065320-turtle.txt
N-Triples: 831-1.0065320-rdf-ntriples.txt
Original Record: 831-1.0065320-source.json
Full Text

Full Text

Shadows Created by anArray of Acoustically Hard CylindersAbove a Rigid PlanebyDavid ClarkB.A.Sc The University of British Columbia, 1992A THESIS SUBMTTr1D IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF ELECTRICAL ENGINEERINGWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAOctober 1994© David Anthony Clark, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.____________________Department of tt1TLLtOi £/L-LrL/L41t/)The University of British ColumbiaVancouver, CanadaDate [7DE-6 (2188)11ABSTRACT“Shadows Created by an Arrayof Acoustically Hard CylindersAbove a Rigid Plane”by David ClarkTransient radiation from, and scattering by, objects of finite size is importantfor applications pertaining to antennas and to target identification. The goalof this thesis work was to study the waves scattered from an array ofcylinders above a rigid plane by running experiments on real models using anacoustic scatter-mapping system developed by Garfield Mellema [12]. Theexperiments were then simulated using MatLab 4.0 on a Sun Sparc 10workstation. The first experiments were run on three different sized cylinderswith a radius of 1.3 cm (0.35 wavelengths ), 5.7 cm (1.54 wavelengths),and 10.6 cm ( 3.15 wavelengths). The second set of experiments were runon a two cylinder array with each cylinder having a radius of 5.7 cm. Theparallel cylinders were placed directly above each other at varying distancesapart. The simulations used the method of images and an algorithm firstpresented by V. Twersky which calculated the incident and scattered fields inthe wave-number domain.In this paper a time domain solution was obtained by applying the inverseFourier transform to the wave-number domain data. Twersky used a modalexpansion to represent the fields scattered by each cylinder. Then theboundary conditions on each cylinder were used to determine each cylinder’s111scattering coefficients. To solve for the boundary conditions on eachscatterer Twersky used the addition theorem for Hankel functions to translatecylindrical co-ordinate systems between cylinders. A problem with thismethod, however, was that it also used a modal expansion which convergedslowly for values of ka much greater than unity (k represents the wavenumber, and a is the cylinder radius). To overcome this convergence problema different translation formula was derived using the integral form of theHankel function and employing the stationary phase approximation. Inaddition, it was found that there was a relation between the scatteringcoefficients of the image cylinders and the scattering coefficients of the realcylindrical scatterers. This result was also substantiated in a study of multiplescattering by two cylinders completed by J. W. Young and J. C. Bertrand [7].Before the field theory is presented, most of the arrival paths are determinedusing appropriate ray paths. Plots of the arrival path length versus speakerposition are presented so that each arrival path on the real data sets and thesimulations can be traced. The results of the real data sets along with theoptical ray path distance plots confirm the correctness of the simulationalgorithm. The simulations correctly predict the arrival distances and thedestructive interference patterns found on the real data sets for the single anddouble cylinder case.ivTABLE OF CONTENTSPageABSTRACT.iiTABLE OF CONTENTS ivLIST OF ILLUSTRATIONS viI. INTRODUCTION 1A. Scattering by Cylindrical Structures 1B. Generation of Real Scattering Data 3C. Thesis Outline 4II. DATA COLLECTION USING AN ACOUSTIC IMAGINGSYSTEM 6A. Description of the Imaging System 6B. Description of the Real Data Sets 10Ill. DETERMINATION OF THE RAY PATH ARRIVALS 11A. Ray Path Arrivals for the Single Cylinder Above a RigidPlane 111. Path Lengths of the Incident and Primary ReflectedArrivals 122. Path Lengths of the Secondary and Tertiary ReflectedArrivals 213. Arrival Distances for the Rays Reflected Twice OffEach Cylinder 26VB. Ray Path Arrivals for the Two Cylinder Array Above aRigid Plane 30IV. FIELD CALCULATIONS AND COMPUTERSIMULATIONS 40A. Derivation of the Fields Scattered by a Cylinder Above aRigid Plane 40B. Derivation of Fields Scattered by Two Cylinders Above aRigid Plane 48C. Description of the Simulation Software 52D. Results for the Single Cylinder Simulation 53E. Results for the Two Cylinder Simulation 55V. RESULTS OF THE REAL DATA SETS 63A. Scattering Results for the Single Cylinder 63B. Scattering Results for Two Cylinder Scatterers 71VI. CONCLUSION 79A. Summary of the Single Cylinder Scattering Results 80B. Summary of the Two Cylinder Array Scattering Results 80APPENDIX A. TRANSLATION DERIVATION FOR THEHANKEL FUNCTION 82APPENDIX B. MATLAB SIMULATION SOFTWARE 89REFERENCES 97viLIST OF ILLUSTRATIONSTables PageTable 3.1 Departure Angles and Corresponding Arrival Angles 35Figures PageFigure 2.1 Block Diagram of the Acoustic Imaging System 7Figure 2.2 Chirp Signal that was Applied to the Speaker 8Figure 2.3 Spectrum of the Chirp Signal 8Figure 2.4 Reference Signal used for Match Filtering 9Figure 2.5 Spectrum of the Reference Signal 9Figure 3.1 Method of Images Applied to the Single Cylinder Problem 12Figure 3.2 Incident Arrivals in the Illuminated and Shadow Regions 13Figure 3.3 Reflection of an Incident Ray by a Cylinder 17Figure 3.4 Secondary Arrivals in the Illuminated and Shadow Regions 22Figure 3.5 Ray Paths of Double Reflection Ray Path Region 27Figure 3.6 Simulated Scattering of the Small Cylinder ( a = 1.3 cm)from Ray Calculations 31Figure 3.7 Simulated Scattering of the Medium Cylinder ( a = 5.7 cm)from Ray Calculations 32Figure 3.8 Simulated Scattering of the Large Cylinder ( a = 10.6 cm)from Ray Calculations 33Figure 3.9 Creeping Waves Departing the Top Cylinder an Arriving onthe Bottom Cylinder 34viiFigure 3.10 Simulated Scattering of Two Cylinders ( d = 1.0 cm) fromRay Calculations 37Figure 3.11 Simulated Scattering of Two Cylinders ( d = 5.0 cm) fromRay Calculations 38Figure 3.12 Simulated Scattering ofTwo Cylinders (d = 25.0 cm)from Ray Calculations 39Figure 4.1 Method of Images Used on an the Single CylinderScattering Problem 41Figure 4.2. Method of Images Used on an Array of Two Cylinders 50Figure 4.3 Simulated Scattering of the Small Cylinder ( a = 1.3 cm)from Field Theory 56Figure 4.4 Simulated Scattering of the Medium Cylinder ( a = 5.7 cm)from Field Theory 57Figure 4.5 Simulated Scattering of the Large Cylinder ( a = 10.6 cm)from Field Theory 58Figure 4.6 Simulated Scattering of Two Cylinders (d = 1.0 cm) fromField Theory 60Figure 4.7 Simulated Scattering of Two Cylinders (d = 5.0 cm) fromField Theory 61Figure 4.8 Simulated Scattering of Two Cylinders (d = 25.0 cm) fromField Theory 62Figure 5.1 Scattering Response of the Small Cylinder ( a = 1.3 cm)from a Real Data Set 65Figure 5.2 Scattering Response of the Medium Cylinder ( a = 5.7 cm)from a Real Data Set 66vuFigure 5.3 Scattering Response of the Large Cylinder ( a = 10.6 cm)from a Real Data Set 67Figure 5.4 Real Signal at x = 125 cm for the Large Cylinder 68Figure 5.5 Simulated Signal at x = 125 cm for the Large Cylinder 68Figure 5.6 Real Signal at x = 95 cm for the Large Cylinder 69Figure 5.7 Simulated Signal at x =95 cm for the Large Cylinder 69Figure 5.8 Real Signal at x = 65 cm for the Large Cylinder 70Figure 5.9 Simulated Signal at x =65 cm for the Large Cylinder 70Figure 5.10 Scattering Response for Two Cylinders (d = 1.0 cm)from a Real Data Set 73Figure 5.11 Scattering Response for Two Cylinders (d = 5.0 cm)from a Real Data Set 74Figure 5.12 Scattering Response for Two Cylinders (d = 25.0 cm)from a Real Data Set 75Figure 5.13 Real Signal at x = 125 cm for Two Cylinders (d = 25 cm ) 76Figure 5.14 Simulated Signal at x = 125 cm for Two Cylinders(d=25cm) 76Figure 5.15 Real Signal at x = 95 cm for Two Cylinders (d = 25 cm ) 77Figure 5.16 Simulated Signal at x = 95 cm for Two Cylinders(d=25cm) 77Figure 5.17 Real Signal at x = 65 cm for Two Cylinders (d = 25 cm ) 78Figure 5.18 Simulated Signal at x = 65 cm for Two Cylinders(d=25cm) 78Figure A.1 Translation of the Hankel Function 83I. INTRODUCTIONScattering by objects of fmite size is important in applications pertaining toantennas, and to target identification. In this thesis the waves scattered bycylindrical structures were investigated using real data sets, field theory, andoptical ray tracing methods to show the existence of creeping waves.Creeping waves are defmed as waves that travel around the circumference ofthe cylinder. Creeping waves scatter in all directions with each scattered raydeparting tangentially from the cylinder. An acoustic scatter-mapping system[121 was used to conduct experiments on six different scattering models. Themodels were then simulated using a time domain solution for the multiplescattering problem. The effects of the rigid plane were included byemploying the method of images.A. Scattering by Cylindrical StructuresAn exact solution for the diffractive scattering of an incident wave by anarray of cylinders was first presented by Twersky in 1952 [1] and [2].Twersky solved the multiple scattering problem by writing the total fieldas an incident field plus various orders of scattered fields. The scatteredfields produced by each cylinder were written as a modal expansion.Twersky then solved for the scattering coefficients by using an iterativeprocedure, which determines the scattering coefficients for m + 1scatterers from the scattering coefficients of the previous m scatterers inthe algorithm. A complete derivation of the iterative procedure can be2found in Chew [3]. A completely general matrix formulation of themultiple scattering equations was first given by Burke, Censor, andTwersky [4] in 1965. V. E. Glazanov [5], R. P. Radlinski and T. J.Meyers [6], and J. W. Young and 3. C. Bertrand [7] also used the matrixformulation method to study the scattered fields produced by a circulararray of cylinders encompassing a single pulsating source cylinder. Inthis thesis the matrix method was used because of the fmdings of J. W.Young and J. C. Bertrand. Young and Bertrand found the matrixinversion method results to be comparable to the iterative approach withthe matrix inversion method being the more efficient algorithm. Studiesthat used the iterative approach include Olaofe [8], in 1970 and Yoursiand Fahy [9] who, in 1973, used this method to calculate the radiationfrom two driven, parallel cylinders. Both the iterative procedure and thematrix inversion method use the addition theorem for Hankel functions.The addition theorem for Hankel functions is the means by which thecylindrical wave functions of one coordinate system can be expressed interms of similar wave functions of a second coordinate system.After the boundary conditions are solved on each cylinder, matrixinversion is used on a truncated set of infmite equations to solve for thescattering coefficients of each cylinder. These scattering coefficients arethen used to calculate the radiation scattered from each cylinder at eachsource and receiver location. The problem with the addition theorem forHankel functions is that the summation converges slowly as the size ofthe cylinder, and/or the frequency contents of the transmitted acoustic3signal increases. Young and Bertrand found that for ka> 1 the additiontheorem for Hankel functions requires N 2ka terms for convergence.To overcome this convergence problem a different translation formula forHankel functions was derived to carry out translations for intermediateand large values of ka. This translation formula was based on a methodpresented by G. N. Watson [101 that was applied to the integralexpression for the Hankel function. The integral expression was thenapproximated using the method of stationary phase shown in [11]. Totransform the data to the time domain, the negative frequency values werecalculated and then the inverse Fourier transform was computed.B. Generation of Real Scattering DataAn acoustic scatter-mapping imaging system designed and built byGarfield Mellema [12] gathered the scattered data and had softwareinstalled so that the scattering results could be plotted in a wiggle tracingformat on a laser printer. The core of the system consisted of a RCElectronics data-acquisition card installed in a XT-compatiblemicrocomputer. The data-acquisition card digitally generated a pulsedsource signal that was sent to a loudspeaker and digitized an analoguesignal received by the microphone. The speaker was mounted on atrolley and its position was controlled by a synchronous motor connectedto the computer via a data-acquisition card. The data was taken at 2 cmsampling intervals over the entire 2.5 m span of the trolley. Data was4collected sixteen times at each sample position to increase the signal-to-noise ratio. After a data set was obtained it was plotted on a laser printer.C. Thesis OutlineChapter I introduces the objective of the thesis and describes how thedata sets were obtained. In addition, chapter I introduces the numericalprocessing used for the simulations and it presents a thesis outline. Inaddition, a discussion is given on previous studies done on multiplescattering by cylindrical structures. Before going into the diffractionproblem presented by curved surfaces in chapters III and W, chapter IIdescribes the acoustic scatter-mapping imaging system used to collect thereal data sets. To create an accurate simulation the physicalcharacteristics of the frequency swept pulse along with other parameterssuch as the sampling rate, number of samples collected, and the type ofwindowing used during fast correlation needed to be known. The fastcorrelation used by the acoustic scatter-mapping system is a method ofcrosscorrelation that multiplies the fast Fourier transformed receivedsignal with a fast Fourier transformed reference signal. Chapter ifidescribes the paths of each arrival that can be seen on the real andsimulated data sets for the one and two cylinder scattering models. Inthis section the optical ray path distance of each arrival is calculated.Plots of each model are included so that each arrival found for thesimulated and real data sets can be traced. In chapter IV the experimentis simulated using field theory in the wave-number domain. The method5of images is used in the derivation to remove the rigid plane, along with adifferent translational form of the Hankel function derived in Appendix A.This translation allows for the scattering coefficients of each cylinder tobe determined using the boundary conditions on each cylinder. After thefields were calculated the data was inverse Fourier transformed to thetime domain. Finally, chapter V presents the real data sets and describesany unpredicted arrivals not seen on the simulated data sets.6II. DATA COLLECTION USING AN ACOUSTIC IMAGINGSYSTEMA study of shadows created by acoustically hard cylinders is demonstratedusing the acoustic imaging lab in the Geophysics building. The acousticscatter-mapping system found in the acoustic imaging lab allowed for thegathering of data from real models to determine how accurate the simulationswere.A. Description of the Acoustic Imaging SystemThe imaging system of Figure 2.1 consists of the following: a loudspeakerthat traverses along a 2.5 m long trolley taking samples every 2 cm via asynchronous motor, a microphone, a D/A converter and an A/D converter,two amplifiers (an audio amplifier to drive the speaker and an electronicamplifier connected to the microphone ), and a 286 XT-minicomputer thatperforms the imaging algorithm and controls the synchronous motor. Thepositioning of the motor is accurate to within 1 mm ( 5 % accuracy). Thetransmitted acoustic signal from the computer to the loudspeaker is a chirpsignal that linearly sweeps over the frequency range from 1-15 kHz. Thecomputer sends a chirp signal every 69 ms at 3 is intervals via a D/Aconverter and an audio amplifier. Figure 2.2 displays the transmitted chirpthat was sent to the speaker and Figure 2.3 shows the spectrum of thechirp signal. The chirp is given a triangular envelope, and has a cosinetaper applied to it in the frequency domain to reduce the sidelobe levels.The dominant wavelength of the transmitted chirp signal is 37 mm. After7the chirp is transmitted the A/D converter takes samples of the receivedsignal (via an electronic amplifier) every 5 jis for 81.92 ms. TheFigure 2.1 Block Diagram of the Acoustic Imaging Systemimaging software then correlates the received signal with the transmittedchirp and truncates the data down to 1860 points so that the data may bestored on a single 1.44 Mb floppy disk. The fast correlation involvestaking a fast Fourier transform of the received signal and multiplying itwith a fast Fourier transformed reference signal. The reference signal is acopy of the received signal that is acquired by directing the speaker into anarea effectively free of scatterers and directing the microphone toward theloudspeaker at a distance of 50 cm. Figure 2.4 displays the receivedreference signal and Figure 2.5 shows the reference signal’sAT. cop.tIb. CopaerI I j:: Sampling Ii ProccwogReceived‘[ CaIculItionZl 11)I.encratio. jAcci signai\8Fig. 2.2 Chirp Signal that was Applied to the Speaker0.01 0.02 0.03 0.04 0.05 0.06F—250020000)0D100050000time(s)1 Fig 2.3 Spectrum of the Applied Chirp Signal430).4.-E1000.075000 10000 15000Frequency ( Hz)9CD 3.5-o0S< 3.U)0a.Sx Fig. 2.4 Reference Signal used for Match FilteringI—2.50 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09time(s)Fig. 2.5 Spectrum of the Reference Signal1510 -5-00 2000 4000 6000 8000 10000frequency ( hz)12000 14000 1600010spectrum. The process repeats itself sixteen times for each speakerposition allowing for averaging of the received data to increase the signal-to-noise ratio.B. Description of the Real Data SetsThe first sets of data were obtained and simulated for a single cylinderabove a rigid plane. The three cylinders that were used as scatterers hada radii of 1.3 cm ( 0.35 wavelengths ), 5.7 cm (1.54 wavelengths ), and10.6 cm ( 3.15 wavelengths). The simplest scattering case, a singlecylinder, confirmed the correctness of the simulations from field theoryby comparison of the results with the real data sets and the ray tracingplots. The next sets of data were collected from a linear array of twocylinders each with a radius of 5.7 cm. The parallel horizontal cylinderswere positioned one above the other. For the three data sets the twocylinders were separated by distances of 1.0 cm ( 0.26 wavelengths),5 cm (1.35 wavelengths ), and 25 cm ( 6.75 wavelengths). For all thesimulations and experiments the transducer was positioned at a constantheight of 23.5 cm from the rigid plane, and 15.5 cm from the center of thecylinder directly above the microphone. The loudspeaker moved alongthe trolley at a constant level of 144 cm from the rigid plane.11Ill. DETERMINATION OF THE RAY PATH ARRIVALSTo provide a means of tracing each of the arrivals and to provide anothermeans of comparison, two MatLab programs were written to compute thelengths of all of the arrivals that could be seen on the simulation and real datawiggle trace plots. From the plots that were generated using these twoprograms, the arrival path of each ray could be detennined.A. Ray Path Arrivals for the Single Cylinder Above a Rigid PlaneThe first program determined the first five sets of arrival lengths for thesingle cylinder above a rigid plane. To convert all of the arrival lengths toarrival times, the lengths were divided by the speed of sound in air at 20°C. Only the first five sets of arrivals were considered because this is thelargest number of visible arrivals on the simulated or real data sets.Depending on the position of the speaker the arrivals came in five sets oftwo or three arrivals per set. Two arrivals came from each side of thecylinder and depending on the speaker position a third arrival came via areflection off the cylinder. For the derivation, the five sets of arrivalswere given the subscripts 1- 5 and the three subsets of the five arrivalswere given the subscripts “a”, “b”, and “c”. The method of images wasused to take into account the reflections off the rigid plane. The methodof images removes the rigid plane and replaces it with a mirror image ofthe scene about the y = 0 plane as shown in Figure 3.1.12o Real Speaker Location y-m X_________VoxReal Cylinder0 Real Microphone LocationLocation of the Rigid Plane hImage Cylinero Image Speaker LocationFigure 3.1 Method of Images Applied to the Single Cylinder ProblemThe problem was split into three by initially considering the first two setsof arrivals that were the result of the incident and primary reflected raypaths.1. Path Lengths of the Incident and Primary Reflected ArrivalsThe incident arrivals are defmed as rays that travel over an obstructedpath, or travel around the cylinder, or are reflected off the cylinder.Primary arrivals are defmed as rays that travel over a path reflectedonce off the rigid plane to the microphone. Fermat’s principle statesthat each ray will travel along a path that makes the ray path distance in13the media stationary. Using this principle the ray paths shown inFigure 3.2 can be constructed. Figure 3.2 shows the incident ray pathsfor the case when the receiver is within the illuminated region of thespeaker and when it is within the shadow region of the speaker. Whenthe microphone is within the shadow region of the loudspeaker thereare two creeping wave arrivals and no incident arrival, and when themicrophone is within the illuminated region of the loudspeaker there isone incident arrival, one creeping wave arrival, and one reflectedarrival.Figure 3.2 Jnëident Arrivals in the Illuminated and theShadow RegionsIlluminated Region for the SpeakerSourceShadow Region for the Speakerarnval_1 bdlamvaU aSourcearrival_i b I jdl /L :1ad2d3 d3Receiverd3Receiver14To find the path lengths of the creeping wave arrivals each arrival pathis sub-divided into the following three parts (i.e., three distancecomponents di, d2, and d3): the incident path from the speaker to apoint tangent to the cylinder, the circular path around the circumferenceof the cylinder, and the departure path to the microphone that alsoleaves a point tangent to the cylinder. Considering only those arrivalsreaching the real receiver point first, the ray path lengths of the firsttwo/three arrivals are determined. A simple substitution is then shownto determine the arrival lengths to the image receiver location. Aspreviously stated, when the microphone is within the shadow region forarrival_la i.e.,I Oi I <atan () where O = atan Yop (3.1)there are two creeping wave arrivals, where 01 is the angle of theloudspeaker with respect to the microphone, a is the radius of thecylinder, x is the current position of the loudspeaker, Yo is the distancefrom the center of the cylinder to the loudspeaker, p< is the distancefrom the center of the cylinder to the microphone, and d3 is thedistance from the microphone to a point tangential to the cylindersurface. One arrival is from the clockwise traveling wave and the otheris from the counter-clockwise traveling wave. To find the path lengthsof these two arrivals the distance from the speaker to the axis of the15cylinder, and its corresponding angle as measured from the cylinder’sx-axis are first determined to be:p = \Jx2+Yo2 (3.2)It X4>=-atan() (3.3)Next the lengths and incident angles of the two ray paths that travel topoints tangential to each side of the cylinder are given by:d1=\[p>2a2 (34)41°a =4> -atan( a=-atan()-atan() (3.5)Ob=4>+atan( a=-ata()+atan() (3.6)where 9a is the incident angle for the clockwise creeping wave, 9b isthe incident angle for the counter-clockwise creeping wave. d1 is thedistance of both rays from the speaker to a point tangential to thecylinder’s edge. Now the departure angles for each creeping wave are16determined, along with the distances of the two departing tangentialrays from the edge of the cylinder to the microphone:d3=\Jp<2a2 (37)O2a+atafl() (3.8)e2b=T-atan() (3.9)°2a is the departure angle for the clockwise traveling wave, e2b is thedeparture angle for the counter-clockwise traveling wave. d3 is thedistance of the departure ray from the cylinder edge to the microphone.The two distances d2a and d2b that the creeping waves travel aroundthe cylinder are given by:d2a=a I °2a°a I (3.10)d2b=aIO2b-ObI (3.11)Thus, the total distance of each incident ray from the speaker to thereal receiver point when the microphone is within the shadow region is:arrival_la = d1 + d2a + d3 (3.12)17arrival_lb = di +d2b + d3 (3.13)When the microphone is located in the illuminated region there is onecreeping wave arrival that can be computed from equation (3.13 ), oneunobstructed arrival whose ray path distance given byarrival_la = /x2+(Yo+p <)2 (3.14)and there is one arrival that is reflected off the cylinder. Figure 3.3shows how an incoming ray will be reflected by the wall of thecylinder. An incoming ray will be reflected normal to a tangent planeto the cylinder wall at the point of incidence.Reflecting CynderFigure 3.3. Reflection of an Incident Ray by a Cylinder18An incoming ray at an angle 0 with respect to the y-axis will arrive atan angle N’ with respect to a normal vector to the cylinder’s edge at thepoint of incidence. The normal vector to the cylinder’s edge is at anangle (p and the outgoing ray is at an angle y from the y-axis. Theincoming ray is perfectly reflected by the cylinder and leaves at anangle N’ with respect to a normal vector to the cylinder’s edge. Fromthe geometry of Figure 3.3 the following equations are obtained:(3.15)(3.16)Substituting equation ( 3.15 ) into equation ( 3.16 ) gives:(3.17)The reflected arrival is first split into distance component 1 and 12. 1 isthe distance of the ray reaching the cylinder from the speaker and 12 isthe distance of the ray reaching the microphone from the cylinder.Equating the vertical and horizontal components of l and 12 gives thefollowing equations:l1cosO=Yo-asinØ (3.18)19l2cosy=p<+asinØ (3.19)l2sin’y=-acosØ (3.20)where is the angle of incidence on the cylinder measured with respectto the cylinder’s x-axis, 0 is the angle of the incident ray, and y is theangle of the outgoing ray. Realizing that the sum of the horizontalcomponents of 1 and 12 must equal the speaker’s horizontal location x,and using equations (3.19 ) and ( 3.20 ) it can be seen that:x =-1 rn 0-12 mn=-(Yo-asinØ)tan0+acosØ (3.21)After substituting equation ( 3.19 ) into equation ( 3.20 ) Newton’smethod is applied to solve for the angles 0 and 0. Newton’s methodprovides a numerical technique to determine the two unknown angles 0and 0. First equations (3.20) and ( 3.21 ) are rewritten as twofunctions of 0 and 0 that sum to zero as shown below:=‘f(0,Ø)=0=acos+(p<+asin)tan’y (3.22)=g(0,)=0=x+(Yo-asinØ)tan0-acos (3.23)20Then the method is iterated through until the two following formulasreach a desired accuracy of 0.5 mm:in±i=øn f(O=On,=4) (3.24)‘gp g’9n+1 =On+ (O=On,Ø=Ø) (3.25)gØ gOwhere 00 and Ø are initial estimates of the variables 0 and theirinitial values are given by:0o=tani(a) (3.26)Po =0 (3.27)where f0, f, gO, and gØ are the differentials of the functions f and gwith respect to the unknown variables 0 and 4. After 0 and 0 aredetermined from Newton’s Method, arrival_ic is given by:21arrival_ic = 11+ 12—Yo-asinØ p< +asinØ— cosO cos( 20-0)(3.28)To obtain the primary reflected arrivals represented as arrivals to theimage receiver point (i.e. the ray paths which reflect off the rigidplane) p < is replaced by Pim in equations (3.1), (3.7), (3.14),(3.19 ), ( 3.22), and ( 3.28 ) to obtain the next set of two/threearrivals, ( arrival_2a, arrival_2b, and arrival_2c ) where pj is shownto be:Pim = p < + 2h, where h = microphone height (3.29)2. Path Lengths of the Secondary and Tertiary Reflected ArrivalsSecondary reflected arrivals are rays which have reflected off the rigidplane and then the cylinder, and tertiary arrivals are rays which havereflected off the cylinder once and the rigid plane twice. First thosearrivals that reach the image receiver point are considered. To obtainthe other two/three arrivals reaching the real microphone point, asimple substitution is again used. To fmd the path of an incident rayreflecting once off the image cylinder Figure 3.4 is considered.22Figure 3.4. Secondary Arrivals in the Illuminated and Shadow RegionsThe arrival that reflects off the real cylinder to the bottom cylinder isincluded in Figure 3.4. Arrival_3c shown on Figure 3.4 was notincluded in the following ray tracing equations or the ray tracing plotsbecause its arrival is not detected on either the real data or simulationwiggle trace plots. Equating the vertical and horizontal components ofl, 12, and 13 the following equations are obtained:Real CylinderIiamval_3aImage Receiver• LocationImage Cylinderl cos 03 = Yo + 2(p <+ h)- a cos 3 (3.30)23l2cosy3=p<-acosØ3 (3.31)l2Sm ‘y3=-asrnØ3 (3.32)where 1 is the length of the incoming ray incident upon the cylinderand 12 is the path length of the outgoing ray that reaches the imagereceiver location. 43 is the angle of incidence on the reflected cylinder,03 is the angle of the incoming ray, and y3 is the angle of the outgoingray. Realizing that the sum of the horizontal components of 1 and 12must equal the speaker’s horizontal location x, and using equations(3.30) and ( 3.32) it can be seen that:x= -1 sin 03 -I- 12 siny3= (Yo+2(p<+h)-acosØ3)tan03+(p<-acos3)tan’y3(3.33)After substituting equation ( 3.31 ) into equation ( 3.32 ) Newton’smethod is applied to solve for the angles 03 and Ø. Newton’s methodgives a numerical method to determine the two unknown angles 03 andØ. First equations (3.32) and ( 3.33 ) are rewritten as two functionsof 03 and Ø that sum to zero as shown below:=f(03,Ø3)=0=asin3+(p<-acosØ3)tany3 (3.34)24= g(8,4) =0 = x + ( Yo + 2(p< + h) -a cos 03 ) tan 03-(p<-acosØ3)tan’y3 (3.35)Then Newton’s method is used until an accuracy of 0.5 mm isobtained. 00 and 4 are initial estimates of the variables 0 and 0 andare given the initial values:Oo=0o=taw1(yO+2+h)) (3.36)After 03 and 4 are determined from Newton’s Method, arrival_3a isgiven by:arrival_3a = l + 12— i Yo+2(p<+h)-acos03 p< -acosØ3— cos 03 + cos( 203 -03)(3.37)Figure 3.4 is considered to obtain the path length of a creeping wavearrival off the top cylinder that reflects once off the bottom cylinderwhen the image microphone is within the shadow region of thespeaker. The vertical component of the creeping wave departure offthe top cylinder is equated to:l1cos 03 =2(p<+h)- a(cos Ø -sin03) (3.38)25Referring to Figure 3.4 and equating the horizontal components of 1and 12 it is found that:a cos 03 = -1 sin 03 +l2srn’y3= (2( p<+h )-a (cos 03-sin 03)) tanO3 + tan (p <-a cos 03)(3.39)After equation ( 3.31 ) is substituted into equation ( 3.32) Newton’smethod is applied to equations (3.32) and ( 3.39 ) Equation ( 3.40) isthen utilized to initialize the values of 03 and Ø as shown on thefollowing page:= =- tanl( 2(p<± h) (3.40)The angle 03 returned by Newton’s method is used to determine thespeaker position at the boundary of the shadow region for arrival_3aand to determine the angle of departure for the two creeping rays offthe cylinder as shown below:03a (3.41)O3btO3 (3.42)2603a is the departure angle for arrival_3a and 03b is the departure anglefor arrival_3b. Newton’s method is used on equation ( 3.32) andequation ( 3.39 ) to initialize the shadow regions before any of the raypath lengths are calculated at the beginning of the two programs thatcomputes all of the path lengths. To obtain the tertiary arrivals definedas ray paths reaching the real receiver point (i.e. arrival_4a) (p - acos ) is replaced by (p +2h - a cos 0) in equations (3.30 ) and(3.33 ) - (3.40) and the exact same operations are carried out to solveforO4, (P4,04a,afldO4b.3. Arrival Lengths for Rays Reflected Twice Off Each CylinderFigure 3.3 is again used to compute the lengths of the two arrivals thatreach the microphone via two reflections off the rigid plane and thecylinder in the illuminated region. For the double reflection 0 is againdefmed as the angle of the incident ray, Øj and 12 are the angles ofincidence on each cylinder, and ‘yj and ‘y are the departure angles offeach cylinder. Referring back to the prior section, any incoming ray atan angle 0 that is incident on a cylinder at an angle 0 results in anoutgoing ray at an angle y ( where y = 20 - 0). Therefore, the outgoingray that is reflected off the bottom cylinder is leaving at an angle yj.To fmd the outgoing angle y Figure 3.3 is used to show that theincoming angle 2 is equal to yj. Making use of equations (3.15 ) -(3.17 ) yields the following relationships:27Figure 3.5. Double Reflection Ray Path= - 02= ‘Yl-02arrival_Sb/.7Real CylinderLocalionarnval5aImage Cylinder=(20i-O1)-02 (3.43)28‘Y2 = 02 - 2jJ=(20i -Oi)-2(20i -Oi -02)=01-2(01-02) (3.44)From Figure 3.5 the double reflection ray path is shown to be made upof three subpaths labeled as l, 12, and 13. Note that arrival_5c shownon Figure 3.5 was not included in the derivation below or the raytracing plots that will be presented later in this section because it wasto weak of an arrival to be detected by either of the real data orsimulation wiggle trace plots. Equating the vertical components ofthese three subpaths gives:l1cosO =(Yo+2(p<+h)-acosOi) (3.45)l2COs•yi =2(p<+h)-a(cos0i +cos02) (3.46)l3cosy=p<-acosØ (3.47)and equating the horizontal components of these subpaths yields:29x=-l1sinO+l2’yi-l3srny2=-(Yo+2(p<+h)-acosOi)tanO+(2(p<+h)-a(cos+cosØ))tany-(p<-acosØ)tany (3.48)a sin Øj = 12 Sill Y1 - 13 sin Y2=(2(p<+h)-a(cosOi +cosO2))tan’yl-(p<-acosØ)tany (3.49)a sin 02= -13 Sm Y2=-(p<-acosØ)tan’y (3.50)To find the three unknowns 9, 0i and 02 Newton’s method is applied toequations (3.48 ) - (3.50). The program that was used to determine arrivalpatterns for the three different sized cylinders created a postscript file foreach size cylinder as shown in Figures 3.6, 3.7 and 3.8. The arrivals for thefirst arrival set are labeled “a”, “b”, and “c” on the plots to indicate whicharrival each hyperbola corresponds to. Arrival “a” is an unobstructed ray or acreeping wave arrival, arrival “b” is a creeping wave, and arrival “c” is a30reflected ray off the wall of the cylinder. On the simulated data sets there aremore than five sets of two arrivals; however, due to reflections off the standof the microphone only five sets of two/three arrivals are seen in the real datawiggle trace plots.B. Ray Path Arrivals for the Two Cylinder Array Above a Rigid PlaneWith the addition of another cylinder comes an additional five sets of twoarrivals; therefore, there are now five sets of four arrivals in total. Sincethe prior section has already dealt with the paths of the rays travelingaround the lower cylinder to the microphone, this section will deal withrays that travel around the top cylinder and arrive at trajectories tangentto the lower cylinder. The insertion of the upper cylinder means thatthree different speaker regions must be considered. Figure 3.9 shows thetwo sets of arrivals and departures for the case when the lower cylinder isin the shadow region. The speaker distance and angle to the top cylinderaxis are shown to be:P>top\JX2+(Yod)2 (3.51)It x>top=atan(yOd) (3.52)where d is the distance between the centers of the top and bottomcylinders. The distance each ray travels to a point tangential to the topcylinder is given by:DistanceTravelled(m)-L00)i\)-LII‘liiIIIIIIIIIIIIIIIIIIIIIr11111111111111111111111111111111111IiIIliiiliiiliii‘IIIliiiIII111liiilIiilIiilIiii,IIIIII0IIIIIIIIIII:(*)IIIIIIIII/IIII??3•)IIIIIIIIIIIIII111111IIIIIiIIIIII\lIIl\1IIIIIIIIII\II\I—4,I0 -‘I-t :r CDCDICl)p)37;--D0o‘<U)//1/13II/‘1i101/1Il?ha1/1JJ11111II1/1/I\—IIIIIIIIIII1II111111II1II111111CDIIIIIIIiIIIIIIIIIIII1111III11IIIIIP0111111111111111111111111111111111111111111111111111111111111/Il111111111111111111111111111111111111111III111111IIIIII1111IIIIIIIIII11111IIIIIIIIIIIIIC”DistanceTravelled(m)-&0UiF\)0101Iii(IiiiilIiiJIII111111111ii11111111111111hhi1:iii1111111111IIIII3iI111111111111111111*IiliiiIiIiiii111111iiiiii111111illIllIiiliiiI\IIliIiIiliiiiiIiIiiiIiIiiIiiiliiiliiiliiiiiiiiiIIIIiIiIiiIiliiiIiiiiIi1iii—IIIiIiljílililtiiiii1i1iiiiCl)liiiiIlllIiIIIiiiIIiilii0iiIII-‘II/’liIiIIII/i1111/I11111Ill/llI)IIiliiiIiI/IIIjliIiI,lii,iIIIIitIii/III/111111!IIII1IiIII(TIiIiliiiifIIilliiliiiliiiIIIIIIi1IliiiiiI0IIliiiIIIIIIII0..IIliiiIIiiliiiIiIIIiII(.1)liiiliii1111/11liiiliiiliiiliii111111111111111111111111/111111/—iiIII1111IIIIIIII1IIIIilI01DistanceTravelled(m)CI\)-(nIllIIIi‘Ii•IIIIIIlII.II crIllliiIIlIIlCiliiiliiililiiiliiiIlliIllIIIIIIIlillIIIill111111IIiiiiIllIIlIIiIII11111IIIl11t1111111IIIIiIiIiIi—‘IiIIIIIIIIIitIIIIIiIIIiiiIIIIIiiiIllIIIillCl)IkiIiIIICiiIIII1 —CI)iIICD-.CDiiii ICCCDCD-‘ o C,) 0 DH1/1IIiiIII//111//’IIIi///3liililtIIII1111III11IIIIIIIlIIIIIIIIliiI’iiI11lI—itiijiiiiiiliiiiiIIIIIIIIIIII1itIIiIi1i1iIIIIIIIIiIIi’i’II0111111it1i11ii1IIIIIIIIiiIIIIIIIII111111111111IIIIiiliii1111111111IIIII1111tilllilttill111111itilIl’illIIIIIIliltII1iil111lIIIICii34(di)top = ‘NIP > top2a (3.53)Figure 3.9. Creeping Waves Departing the Top Cylinder and Arriving onthe Bottom Cylinderand the arrival angles of each ray are found to be:x (di)t0(0a)top = - atan ( Yo-d - atan a ‘Top CylinderFour DepartingCreeping Waves off theTop CylinderBottom Cylinder(3.54)35(0b)top = - atan Yo-d + atan ((di)top) (355)When the speaker is completely within the shadow region there are fourcreeping wave departures off the top cylinder that arrive tangential to thelower cylinder. The departure angles and corresponding arrival angles foreach possible ray path are given in Table 3.1 as shown below:Table 3.1. Departure Angles and Corresponding Arrival AnglesDeparture Angle off the Top Arrival Angle on the BottomCylinder Cylinder0 0It12a 12a-snr’()2a 2asur1()Once the departure angles for the upper cylinder and the arrival angles forthe lower cylinder are determined, the method shown in section llI.A. isrepeated to find the lengths of the of the five possible paths for eachcreeping wave on the lower cylinder. Figures 3.10, 3.11, and 3.12 showthe results of the two cylinder array scattering problem. The arrivals inthe first arrival set are again labeled from “a” to “e” to correspond withone arrival path: arrival “a” is an unobstructed ray, a creeping wave, or36a double creeping wave that travels vertically between the two cylinders,arrival “b” is a creeping wave or a double creeping wave that travelsvertically between the two cylinders, arrival “c” is a reflected ray off thebottom cylinder, arrival “d” is a double creeping wave that travels eithervertically or diagonally between the two cylinders, and arrival “e” is adouble creeping wave that travels diagonally between the two cylinders.DistanceTravelled(m)r\)-&Q•C)01I\)-01IIIIIIliiiliiiIIIIIIIIIIIliiiIIilillIyiIiiIIIIIiI:IIir’IIIIIIIIII111111:‘11111111:11,liiiI(JI11111IllIllIllIllIllCIj1.IIIIl—C)IIliiiliiiIiliiiIIIIIIIII111111IIliii111111IIliiiIiIiIiIi1IiCIIIIIIIIIIIII111111111IIIl;IIiIIII:i11111IIIiIII\111:I\\IIIIii1IIIilIihIhIICl)IliIIhI—1HH11111111hII/Ii1)iihI/hi/1oiihhhhIiiIII11111I1111113/:1iuu/hIIi,‘hih1hlI?hiIIh1111h1i1111111lbIIIIibilliIil11111111IIIIIIIIIIliiiliiiiii11i11iiiIIhihiI1iI1II1iiIIIII1h1IilhillIIIIIruiIhIIIIhh11Ilhh11111Iii?LIhhhIhI11h11111111ii1jiilliiiIilIIIiiiIIiiiIIIiiII—DistanceTravelled(m)0_C)f’)IIIIIII1IliiiIIIIIliii:111111:1II1111::,IIiIiIIiIIIIIIII11111IIIIliiiI111111111111‘Iil11111IIIIIIIlIiIIiI?IiIIIII11111III\IIIII(71‘IIIIIIIIIIiIiIII1111IIIIIiIIIIIii—..IIII,IiIIIIII\III11111IiiiiIIIIIIIIIII1IfiIi\Iii—.IIIIIIIIIIliiiIlIIfII(1IIIIIIIIIIII0CI)I\II\IIii11 ICDII \IIII0i’/IHI0III0\IC)5(I)/III1//II’/0Dlii31////-IIII1IfIl//// /I)/i‘‘ili/—..‘III•1II•f•ijIIIIfIJI /11III0i/1fflfiff1ff13I/II///1/11lii/Il1IIIIIII11/iiI’fI1fI1II/ffI1f1/1fIIII//IJ/I/fII/I1111111CDI-/11’II111II1IIIIII:i’II:f1IIIIIII:IIIII0!‘L11’I1f111If1IIII111f1111111C),00DistanceTravelled(m)F\)-0(A)(fl1\)(31III111111111111IiIiiIIllII::lI,11111111IIllt1.1h:hi,1b.lIIIIlIIIIIIIllIIllIIIIIIIIIIIIII‘IIIII—IiIiII:CA)IIIliltiIilIiiiiiiiiIiiIiiliiiIiiiiIIIliiIllillIIlIlilillillI1l1l1lliii1iiikIIiiiiIIII\IIIliiiIIliltIiiIiIi111111lilIilkllilllIiliiii1111111—.1..liiiIiliiiIiliiiiiliiiII1111101111111IIliiiIIliiiIllV¶IliiiIIIilk0)11I0CDII\it(_)—‘¶IIII—Jk’iioI)lI11(1)(tli‘III Ii IiII¶/fD111/)-.IIiI1Q3—ii11111lI¶1111iIi—C)0“‘I/hiiiliiilf31111iiIlIIiii1111I11111liiIIllIIlI1IiIlIIii1itiii1I10II/IlIiiJiitiiiiillIIIII/11il111113II/i1i1¶1iiI’II1li/iliIliltlthihii/111111110)iiIiI1IIII11/1II11/I1I111/iIi1ilil1IIii1(I.IIJ)i/‘1iiIIII1iilIIIIIIilitiiiiiiiIIIIIliiiIIIiliiiIIliiiIIIIIiIIII0IiIIii11111tIIlIlIIIIIItI?IiiliII1iliiiiiIiC))IIIIhI11111liiiII11111IllliiiliiiIIIlIltI•liIliltl1\)Ill111111llIfl111111111111III(3140IV. FIELD CALCULATIONS AND COMPUTER SIMULATIONSTwo computer simulations were written using MatLab 4.0. The first programwritten was a simulation of the single cylinder scatterer above a rigid plane,and the second program simulated the scattered arrivals for a two cylinderarray above a rigid plane. Appendix B gives a listing of the MatLab sourcecode for the single cylinder scattering simulation along with the MatLabsource code for the Hankel function translation. The simulations tookanywhere from 3-10 hours to run on a Sparc 10 workstation with the amountof time increasing as the size, or number of cylinders increased.A. Derivation of Field Scattered by a Cylinder above Rigid PlaneIt was desired to determine the field created by an incident line sourceand a single cylindrical scatterer above a rigid plane at a point that isdirectly beneath the cylinder. Using the method of images the rigid planewas removed and a mirror image of the problem was created about they = 0 plane as shown again in Figure 4.1. The image source was equal tothe real source in magnitude and phase so that the derivative of the fieldwent to zero along the rigid plane. The one cylinder problem was thentransformed into a two cylinder scattering problem with two sourcelocations. A method presented by Burke, Censor, and Twersky, andillustrated by Constantine Balanis [13] was used to generate an equation41Figure 4.1. Method of Images Used on an the Single Cylinder ScatteringProblemfor the fields scattered by the two cylinders. The method starts by notingthat the incident field behaves as a magnetic line source (TMpolarization) because the normal derivative of the field goes to zero at thesurface of the cylinder. The incident field of each source is represented asa line source as shown below:P1 = Ho(2) (k I p-p’ I )Real Sourcep1.Real Cylinderp1x1 rnReceiverp2Image CylinderImage Source(4.1)42where p is the distance from the center of one of the cylinders to thetransducer, p’ is the distance from the center of the cylinder in question tothe speaker, k is the wave number, and Ho(2) represents a zero orderHankel function of the second kind. Using the addition theorem forHankel functions taken from [14], the incident field due to the real andimage sources is re-written as a modal expansion as shown below:P = Jn(kp1)Hn(2)(kp1texp( -jn(a1a’)) (4.2)1im Jn(kp1)Hn(2)(kpim’) exp( -jn(a1cx’)) (4.3)The above equation is written in a cylindrical coordinate system with itsaxis centered at cylinder 1. Pi P1t, and P’im are the distances from theaxis of cylinder 1 to the microphone, real speaker, and the image speaker.a1, x1’ afld aim’ are the corresponding angles of the microphone, realspeaker, and image speaker with respect to the axis of cylinder 1. Thescattered fields created by each cylinder at the receiver point M are alsowritten as modal expansions given by:P1 = An1 Hn(2)(kp1exp( -jna1) (4.4)= An2Hn(2)(kp exp( -jna2) (4.5)43where the An’s are used to represent the scattering coefficients of eachcylinder. Thus the total field at the receiver point M can be written bysumming equations (4.2), (4.3 ), (4.4 ), and (4.5 ) as shown below:1m = Jn(kp1)Hn(2)(kp1‘) exp( -jn(a1-a1’))+ Jn(kp1)Hn(2)(kpim’) exp( -jn(a1 -aim’))+ An1 Hn(2)(kp1exp( -jna1)+ An2Hn(2)(kp exp( -jna2) (4.6)To find the scattering coefficients the boundary condition which statesthat the normal derivative of the field with respect to p goes to zero oneach cylinder is applied. Using this boundary condition on cylinder 1results in:Ti=a o (4.7)Before the boundary condition of equation (4.7 ) can be applied thefields of cylinder 2 need to be written in a cylindrical coordinate systemwith its axis centered at the axis of cylinder 1. This is accomplished by44applying a translation theorem for Hankel functions that is presented inAppendix A. This method the derivation of a formula for the derivativeof a translated Hankel function using the stationary phase evaluation of anintegral and a method presented by G. N. Watson. The derivation startsby using the Hankel function’s integral representation. After translatingthe fields of cylinder 2 to a cylindrical coordinate system centered atcylinder 1 equation (4.6) is differentiated with respect to Pi• Then Pi isset equal to a (where a = cylinder 1 radius ), and then the expression forthe total field is set equal to zero as shown below:Ip=a =0 0 Jn’(ka) Hn(2)(kp1exp( -jn(a1a))+ Jn’(ka) Hn(2)(kpim’) exp( -jn(cxi-cxim’))+ An’ Hn(2)’(ka) exp( -jna1)+ An2 (Hn(2)(kpexp -ma2)) Ip=a(4.8)Equation (4.8 ) gives 2n + 1 linear equations in terms of the uiiknown1 2coefficients An and A . To solve for both coefficients another 2n + 1linear equations are required. The other 2n + 1 linear equations comefrom the boundary condition which states that the derivative of the fieldon the second cylinder must also go to zero. In order to apply this45boundary condition the translation method shown in Appendix A is usedagain to write the derivative of the scattered fields produced by cylinder 1in a cylindrical coordinate system with its axis centered at the axis ofcylinder 2. Taking the derivative of equation (4.6 ) with respect to P2setting P2 equal to the radius of the second cylinder ( P2 = a), andequating the entire expression to zero gives the following expression:ap2 1p2=a = 0zz 0 = Jn’(ka) Hn(2)(kptexp( -jn(a2a’))+ Jn’(ka) Hn(2)(kPm’) exp( -jfl(a2a’))+ A’ (Hn(2)(kp1exp -jna1))1p2=a+ An2Hn(2)’(ka) exp( -jncz2) (4.9)To solve the linear system defined by equations (4.8 ) and ( 4.9) amatrix equation is created to solve for the unknown scatteringcoefficients. This is accomplished by rewriting equations (4.8 ) and(4.9 in a different form as shown:ömnAn1Hn(2)’(ka) exp( -jrnx1)+ Lmn’ An2 = - iil (4.10)6mn An2Hn(2)’(ka) exp( -jna2)+ Lmn2An’ = - (4.11)46where ö is the unit delta operator (i.e., & = 1 when m = n), and thevariables 141n1 and lyn2 are defined by:Jn’(ka) H(kp1exp(-jnx1)+ Jn’(ka) Hn(kpim’) exp(-jna1’)(4.12)Jn’(ka) FL(kp2texp(-jna2’)+ Jn’(ka) H(kp21exp(-jna21)(4.13)and the variables Lm1 and L2 are defined by:(Hn(2)(kpexp -jnc)) Ip=a (4.14)(Hn(2)(kp1exp -jna1)) p2a (4.15)Finally using equations (4.10) and ( 4.11) gives the following matrixequation that is used to determine the unknown scattering coefficients.Rewriting equations (4.10) and ( 4.11) in matrix form yields:I Hn(2)’(ka) exp( -jna1) L’ II A1 I = I’4JnhlI Lmn2 Hn(2Y(ka) exp( -ma2) i A112 I I2I(4.16)47Therefore,I An1 I = I Hn(2)’(ka) exp( -jna1) L1 i-i I NJn1 II An2 I I Lmn2 Hn(2)’(ka) exp( -jna2) I I-ijJn2’I(4.17)Once all of the unknown scattering coefficients are determined usingequation ( 4.17 ) they are substituted back into equation (4.6 ) todetermine the total field at the desired receiver location M. Because thespeaker used to obtain the real data sets behaves like a point source anadditional multiplication term is included in each of the summations inequation ( 4.6). This scaling term, an approximate correction, is foundusing the asymptotic formula for the Hankel function along with themathematical representation of a point source. Because the speaker isalways in the far field the line source can be represented as:lme = (4.18)The mathematical representation for a point source is shown below:—jkpP0j= e (4.19)4itpTo approximately convert the line source to a point source equation(4.18 ) needs to be multiplied by the following scaling factor:481 4 I k -iPsce —-I—-—e (4.20)j v 32itEquation (4.20) is first applied to equation (4.17 ) without the term.Then after the simulated field data is inverse Fourier transformed theterm is included by realizing that for each data position p = c n Atwhere c is the speed of sound, n is the sample number and At is thesampling interval.B. Derivation of Fields Scattered by Two Cylinders Above a RigidPlaneThe derivation for the total field at any point M for the two cylinder arrayis similar to the derivation for the single cylinder except that now thereare two more cylinders for which scattering coefficients need to bedetermined. First the method of images is used to remove the rigid planeand create a mirror image of the scatterers and the speaker below therigid plane as shown on the following page in Figure 4.2. Because of thenumber of translations that are necessary when using Twersky’salgorithm, it was decided to run the next experiment on a two cylinderarray with one cylinder being placed directly above the other. Theexperiments were then completed for various spacings between thecenters of the two cylinders. The sources are again assumed to be line49sources and then the scaling term described by equation ( 4.20) wasapplied to convert the sources to a point source. To find an expressionfor the total field the contributions of the incident field plus the scatteredfields are summed together as shown below:Jn(kp1)Hn(ktexp jn(-)++ AlH(kp1)exp jna+ A2H(kp)exp(jna2)+ A3H(kp)exp(jna+ An4 H( kp4 ) exp(jn4 (4.21)50Figure 4.2. Method of Images Used on an Array of Two Cylinderswhere Pi P2 p3 and p4 are the distance from each cylinders center to themicrophone, and (X1, a2, (X3 and a4 are their corresponding angles.Using the method shown in Appendix A the derivative of the total field iswritten in a cylindrical coordinate system centered at the axis of thecylinder in question. Then the boundary condition which states that thenormal derivative of the field on each of the cylinders goes to zero isapplied. After equating the derivative of the total field to zero a solutionto one of the scattering coefficients in terms of the other three scal±eringcoefficients is obtained. The process of solving for the boundaryp1Sourceal’p1Top CylinderBottom Cylinderp2Real ReceiverImage BottomImage Top Cylinderp3Image Source51conditions on the other two cylinders is repeated to obtain 3(2n+1) linearequations in terms of the unknown scattering coefficients as given below:IHn(2)’(ka)ei1j1 L3I him4 Hn(2Y(ka)j2 Linn5 him6I Lm2 Lnn8 Hn(2)’(ka) e-J’3 him9I jpj12 Hn(Y(ka) -jncx4iiA1 I I_141n11i &2 i =IlA3 I-NJn3iIlA4 I-N1n4l(4.22)where the L and qJn are defined in a similar way as for the two cylindercase. The scattering coefficient vector is now found from the followingmatrix equation:I An’ I I Hn(2)’(ka) jn(X1 J1 i2 i3i I = I h4 Hn(2)’(ka) jna2 j5I An3 I I him7 Lmn8 Hn(2)’(ka) ein3 j9I An4 I hmjO [1 1 j12 Hn(2Y(ka) e-nO4(4.23)Once the scattering coefficients have been determined they are insertedback into (4.21) to obtain the total field at the desired receiver locationM.I’I-NJnhII Iin2II I-’qifl3I I-N1n4152C. Description of the Simulation SoftwareAt the beginning of each simulation an initialization routine was run topredetermine all of the Bessel and Hankel functions that were used, alongwith the distances from each cylinder’s central axis to the microphone’sposition and their corresponding angles. The distances from eachcylinder’s axis to the speaker’s position and their corresponding angles,along with their associated Hankel functions were calculated in the outermost ioop which iterated over the 125 speaker positions. The inner loopcovered the frequency range covered by the chirp (i.e. from 1000-15000Hz). The simulations used a sampling frequency of 102.4 kHz and theiterations in the frequency domain were done every 50 Hz102.4kHz2048 points ). This frequency samplmg mterval gave 281 frequencypoints over the frequency band of the speaker and 6.62 m of data. Thesimulated data used 2048 points at each speaker position to cover thefrequency band; however, only 1024 points were saved. The 1024 pointsthat were discarded contained the weak arrivals which had travelledfarther than 3.31 m. The simulations have a resolution of 3.4 mm and thereal data sets have a resolution of 1.7 mm. The modal expansions carriedout in the simulations were found to converge for orders of n equal toka+ 1. Before taking the inverse Fourier transform, the negativefrequencies were filled in by noting that for real data, the values at thenegative frequencies are the complex conjugates of the values for thepositive frequencies. In addition, the data was multiplied by thespectrum of the acmal transmitted chirp and match filtered using the same53reference signal that was used to obtain the real data. After filling in thefrequencies from -Fs to Fs (where Fs equals the sampling frequency ) thedata was inverse Fourier transformed to give a 2-dimensional position vs.arrival time/distance plot. To make the comparison of the simulated andreal data sets easier, the simulated data sets were plotted in the samewiggle trace format as the real data sets. This was achieved byconverting the floating point MatLab data to 16 bit unsigned integer sothat the plots of the simulation results could be transformed into apostscript wiggle trace using the software provided by Garfield Mellema[12]. The correctness of the simulation for a single cylindrical scattererwas checked by using some previously plotted results in Bowman andSenior [15]. This was achieved by moving the single cylinder away fromthe rigid plane so that the effects of the rigid plane became negligible.The results plotted in Bowman and Senior displayed the scattering patternof a single cylinder as a function of the receiver angle for varying valuesof ka, with the source located on the y-axis directly above the cylinder.D. Results for the Single Cylinder SimulationThe results of the simulations for the single cylinder case are shown inFigures 4.3, 4.4, and 4.5. The arrivals in the first arrival set are labeled“a”, “b”, and “c” so a comparison can be made with the ray tracingsimulations. From chapter 3, arrival “a” is an unobstructed or creepingwave arrival, arrival “b” is a creeping wave arrival, and arrival “c” is areflected ray off the cylinder’s wall. The simulation results from field54theory ccincided with the arrivals predicted from the optical ray pathmethod of the prior section. As the cylinder size increased the number ofvisible arrivals increased because of the increased number of reflectionsoff the cylinder’s wall. The smallest cylinder plot shows four sets oftwo/three arrivals, the medium cylinder plot has six sets of two/threearrivals, and the large cylinder plot has seven sets of two/three arrivals.In addition, as the cylinder radius increases the arrivals due to creepingwaves become more visible because the two arrivals in each set traveldistances around the cylinder that differ more in length when the speakeris at the ends of the trolley. The arrivals due to reflections off thecylinder’s wall are best seen in the medium cylinder plot for arrival settwo and the small cylinder plot for arrival set one. Destructiveinterference patterns can also be seen on the scattering plots for themedium and large cylinder. Destructive interference occurs for bothcylinders at two locations close to where the speaker is at its closestrange. These destructive interference patterns occur for every arrivalseen on the wiggle trace plots when the speaker is within this region. Thedestructive interference patterns are caused by the incident wavesdestructively interfering with the creeping waves coming off the oppositeside of the cylinder. Destructive interference occurs when the path lengthdifference between the incident waves and the creeping waves is equal to- ( 18.5n mm in distance where n = 1, 2,...). After the simulationswere complete it was found that the scattering coefficients of the secondimage cylinder could be determined from the scattering coefficients of thereal cylinder by utilizing following relation:55(4.24)This result was also shown in a previous study of the multiple scatteringby two cylinders carried out by Young and Bertrand and is a result of thesymmetry of the problem after the method of images is applied. Equation(4.24) shows how the scattering coefficients of the image cylinders canbe determined from the real cylinder scatterers. Making use of thisrelation would reduce the number of equations (and more importantly thenumber of translations ) needed to solve for the scattering coefficients bya factor of two. In addition, it would decrease the amount computationtime considerably because of the large amount of floating pointoperations required to compute the inverse of matrices.E. Results for the Two Cylinder Scatterer SimulationThe simulation results for the two cylinder array problem are given inFigures 4.6, 4.7, and 4.8. The five arrivals in the first arrival set arelabelled from “a” to “e” so that the labelling on the simulations match thelabelling on the ray tracing plots. To determine which arrival each labelrefers to see chapter 3, section B. The simulation results for the twocylinder array show the same destructive interference patterns as for thesingle cylinder. As the distance between the centers of the cylindersincreases, the size of the destructive interference region decreases. Whenthe cylinders are close together four arrivals can be seen in the firstc1WH1tBillTflT1(.4 (.8 0 00-(.4oio2030405060708090100110120130140150160170180190200210220230240C’SpeakerPosition(cm)Figure4.3SimulatedScatteringbytheSmallCylinder (a—1.3cm)fromFieldTheorya00-0 00-8-IIIIIIIIIIIIIIIIIIoio2030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure4.4SimulatedScatteringbytheMediumCylinder (a=5.7cm)fromFieldTheory0 0 00 0 ClIIIIIIBoio2030405060708090100110120130140150160170180190200210220230240Speaker Position(cm)Figure4.5SimulatedScatteringbytheLargeCylinder (a=lO.6cm) fromFieldTheoiy59four/five sets of arrivals depending upon which location in the plot isbeing investigated. Near the edges of the plots only four of the possiblefive arrivals are seen in each arrival set because only one creeping wavefrom the top cylinder is detected by the microphone. Closer to the middlespeaker location four arrivals can be seen in each arrival set due to thetwo creeping waves coming from the top cylinder. The third and fourthdetected arrivals are the result of creeping waves coming from the topcylinder traveling along a path from the top cylinder directly downward tothe bottom cylinder. The creeping waves leaving the top cylinder thattravel diagonally between the two cylinders are not detected because theymust travel a farther distance around the top cylinder, thus becomingmuch more attenuated. The arrival vs. speaker location plot for the casewhen the cylinders are separated by 5.0 cm shows the first four arrivalsmore clearly because the arrivals have become more separated. As thedistance between the cylinders is increased to 5.0 cm and 25.0 cm arrival“d” becomes more separated from the preceding three arrivals at the sidesof the plots because of the larger change in ray path distance with respectto speaker location. In addition, arrival “e” becomes less separated fromthe preceding arrivals due to a smaller change in ray path distance for allof the speaker positions.C C 00111U11UWThBWfl1WaC—dbC.’)-0102030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure4.6SimulatedScatteringbyTwoCylinders(d=1.0cm)fromFieldTheoryC 1 C-00C.)—4—e_______________IC00-C)- C 00-—0102030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure4.7SimulatedScatteringbyTwoCylinders(d=5.Ocm)fromFieldTheoryC 004(IiIII!I1rI(l)LJ.1)111111)l_.J141144( IillIII 41,..‘113}fllfdlL’1ItI1U1flJiiJi.11_L1.____________C.)iii,(i1- C 00-—0102030405060708090100110120130140150160170180190200210220230240Speaker Position(cm)Figure4.8SimulatedScatteringbyTwoCylinders(d=25cm)fromFieldTheory63V. RESULTS OF THE REAL DATA SETSThe six real data sets all agreed with the predicted arrivals of the ray pathmethod and the field calculations method; however, there are unpredictedarrivals on the real data sets. These arrivals are the result of the microphonenot being a perfect point target, and are also the result of the experiment notbeing done in a room with non-reflective walls. The results of themicrophone not being a perfect point target are repeated arrivals instead ofone single well defined arrival. Because of the microphone’s physical sizethere is a “ringing” effect around each arrival because the main arrival isdetected in the middle of many other arrivals. The arrivals caused by thewalls are not as detrimental to the results because they show up as a weak setof two diagonal lines near the bottom corners of the plots. All of the data setsare labelled in a similar way as the simulation and ray tracing plots to allowsome insight into each arrival’s path.A. Scattering Results For the Single CylinderThe wiggle trace plots from the real data sets for the single cylinderscatterer are given in Figures 5.1, 5.2, and 5.3. The plots have their firstthree arrivals labelled “a”, “b”, and “c” just as for the ray tracing plotsand the simulation wiggle trace plots. In addition, Figures 5.4 - 5.9 showsingle traces for the real and simulated large cylinder data at threedifferent speaker positions. The single trace plots have the first three64arrival sets labelled as “1”, “2”, and “3”. Arrival set “1” is theunobstructed and creeping wave arrivals and the reflected arrival off thecylinder. The main arrival sets in the single trace plots for the largecylinder scatterer all have the same time delay. A point of interestthough, is that the amplitude of the second arrival set is larger than theamplitude of the third arrival set even though this set has travelled afurther distance. Arrival set “2” is the rigid plane reflection arrivals, andarrival set “3” is the double reflection arrivals off the rigid plane and thenthe cylinder. The scattering results from the real data sets show three setsof two/three arrivals for the small cylinder, four sets of two/three arrivalsfor the medium cylinder, and six sets of two/three arrivals for the largecylinder. For all three cylinder sizes, the number of arrivals detected wasless than the number of arrivals found on the three simulations. Thereason for the decreased number of arrivals was that each cylinder’s wallwas not a perfect acoustically hard surface. Some of the incident energywas transmitted through the cylinder, thus not all of the energy wasreflected. The real data sets for the medium and large cylinder displaythe shadow regions created by destructive interference. The datacollected for the medium cylinder shows two shadow regions and the datacollected for the large cylinder shows four shadow regions. Two of theseshadow regions occur close to the middle speaker location and occur forevery arrival seen on the real data sets. The extra two destructiveinterference shadows for the large cylinder occur for the creeping wavearrival in the first set of two arrivals.0 (.4 0 00 00-oio2030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure5.1ScatteringResponseoftheSmallCylinder (a-4.3cm)fromaRealDataSet0 0= E-o 00-oio2030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure5.2ScatteringResponseof theMediumCylinder (a—5.7cm)fromaReal DataSet0 (.4-0 0011JJ0‘CC00-8—(.4 (.4 0 00-(.4—0102030405060708090100110120130140150160170180190200210220230240Speaker Position(cm)Figure5.3ScatteringResponseof theLargeCylinder (a=1O.6cm) fromaReal DataSet6800 0.00115x10°1-Fig. 5.4 Real Signal at x=125 cm for the Large Cyl.86a)4-.0.E20.002 0.003 0.004 0.005 0.006 0.007 0.008time(s)Fig. 5.5 Simulated Signal at x=125 cm for the Large Cyl.0.009 0.01I I I I I I I I•04-0.E 0—0.5—10I I I I I0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009time (s)0.0169U)V0.E—10 0.001 0.002 0.003 0.004 0.005time (s)0.006 0.007 0.008 0.009 0.011 Fig. 5.6 Real Signal at x=95 cm for the Large Cyl.0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01time (s)1010 Fig. 5.7 Simulated Signal at x=95 cm for the Large Cyl.a)V4-0.E1.5 I1-0.50w-Jv—0.5-I I I IU)0.E1 Fig. 5.8 Real Signal at x=65 cm for the Large Cyl.87060).4-’0.E2ftvvI I I II I I I0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009time(s)1 010 Fig. 5.9 Simulated Signal at x=65 cm for the Large Cyl.0.01I .ij10.50—0.5—1I I I II‘v4’÷0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009time(s)0.0171B. Scattering Results for Two Cylinder ScatterersThe wiggle trace plots from the real data sets for the two cylinder arrayexperiment are given in Figures 5.10, 5.11, and 5.12. The real data plotsare have three of their first five arrivals labelled from “a” to “c”. Thelabelling was done in the same fashion as for the two cylinder simulationplots and the two cylinder ray tracing plots. Arrival “a” corresponds toan unobstructed ray, a creeping wave, or a double creeping wavetravelling vertically between the two cylinders. Arrival “b” is a creepingwave or a double creeping wave that travels vertically between the twocylinders. Arrival “c” is the reflected ray off the bottom cylinder. Onlythree of the five arrivals are labelled because the other two arrivalscoming from the top cylinder are undetectable from the real data plotspresented. Figures 5.13 - 5.18 show single trace plots for the real andsimulated two cylinder data for various speaker positions. The first threearrival sets are labelled “1”, “2”, and “3” where arrival set “1” is theincident arrivals, arrival set “2” is the reflected arrival set off the rigidplane, and arrival set “3” is the double reflected arrival set off the rigidplane and then off the bottom cylinder. The time delays for the first threearrival sets on the real data single trace plots all agreed with thesimulation single trace plots; however some of the later arrival sets afterarrival set “3” have different time delays. One notable differencebetween the single trace plots for the real and simulated data is that thesecond arrival set for the real data is much lower in amplitude than the72simulations second arrival set. The scattering results for the two cylinderarray from the real data sets did not detect any of the creeping wavearrivals coming off the top cylinder. The “ringing” that occurred aroundeach arrival set was much too large in intensity to detect arrivals “d” and“e” in any of the wiggle trace plots. To detect these arrivals a smallermicrophone would need to be used to reduce the “ringing” effect. Inaddition, the destructive interference patterns did not change in size as thedistance between the cylinders was changed. The arrivals seen onFigures 5.10,5.11, and 5.12 are almost identical to those arrivals seen onthe single medium cylinder wiggle trace plot ( Figure 5.2). Thus, theaddition of the top cylinder did not add any new detected arrivals.C C C 00 •tI C — 0 102030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure5.10ScatteringResponseofTwoCylinders(d1.0cm) FromaReal DataSet00,-‘0 a.)00-c4-0102030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure5.11ScatteringResponseof TwoCylinders(d—5.Ocm)FromaRealDataSetC C-00 C‘C C C—0102030405060708090100110120130140150160170180190200210220230240SpeakerPosition(cm)Figure5.12ScatteringResponsefor TwoCylinders(d=25cm)FromaRealDataSet76a,4-05X Fig. 5.13 Real Signal at x=125 cm for Two Cyls. (d=25 cm)4-a)V4- -0.E2-I I I1-0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01time(s)1010 Fig. 5.14 Simulated Signal at x=125 cm for Two Cyls. (d=25 cm)1.5 I I I I2’1-0.5 -—0.5 -—1 II I I0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01time(s)77a)-D=4-a.EFig. 5.15 Real Signal at x=95 cm for Two Cyls. (d=25 cm)54a)4-a.E210 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01time(s)1 010 Fig. 5.16 Simulated Signal at x=95 cm for Two Cyls. (d=25 cm)21-O—JN—1 —012I I I I I I I I I0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009time(s)0.0178a)V4-Ei Fig. 5.17 Real Signal at x=65 cm for Two Cyls. (d=25 cm)65a)V4-0.E3210time(s)Fig. 5.18 Simulated Signal at x=65 cm for Two Cyls. (d=25 cm)0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01I I I I I I21-0—1 —01 2I I I I I I I0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009time(s)0.0179VI. CONCLUSIONThe simulation algorithm using the matrix inversion method derived byBurke, Censor, and Twersky and employed by Glazanov and Radlinski andMeyers gives a method of obtaining the arrival patterns from two-dimensionalscattering models exhibiting circular symmetry. The required translationsused to solve for the scattering coefficients can be carried out using theaddition theorem for Hankel functions at low frequencies or the methodpresented in Appendix A for intermediate to high frequencies. Models whichinclude a rigid plane can use the method of images to remove the rigid planefrom the problem by replacing it with a mirror image of the scene on the otherside of the plane. After the method of images is employed the scatteringcoefficients of the image cylinders can be found from the scatteringcoefficients of the real cylinders using the relation presented in chapter 4.Finally, the calculated fields can be mapped from the wave number domain tothe time domain by using symmetry arguments in the wave number domainand then calculating the inverse Fourier transform. Two sources of error inthe field theory simulations are the radiation pattern of the speaker and theradiation pattern of the microphone were not included in the simulations.Sources of error in the real data plots are the spurious arrivals from objectsother than the walls of the building.80A. Summary of the Single Cylinder Scattering ResultsComparison of the simulated single cylinder wiggle trace plots with thereal data single cylinder wiggle trace plots proves the veracity of thesimulation algorithm. The locations and intensities of the hyperbolapatterns found on each of the three simulations for the single cylinderscattering problem all agreed with those arrivals found on the real datasets. The differences between the simulations and the real data sets werethe extra arrival sets found on the simulation plots and the wallreflections. The reason for the decreased number of arrivals on the realdata was the spurious multiple reflections off other objects within the lab.Finally, the locations of the hyperbolas all agreed with the ray tracingplots of chapter 3 which used optical ray paths to calculate the distancesof most of the possible ray paths.B. Summary of the Two Cylinder Array Scattering ResultsThe results for the two cylinder scattering problem from the simulationsshowed that the creeping waves travelling around the top cylinder andthen directly downwards to the bottom cylinder would be detected.Creeping waves traveling diagonally between the two cylinders were tooattenuated to be detected by the microphone, or the simulations. The twocylinder simulations also showed that as the distance between the twocylinders increased, the size of destructive interference regions foundnear the center speaker location decreased. The wiggle trace plots from81the real data sets for the array of two cylindrical scatterers did not showany change at all as the distance between the cylinders was changed. Inaddition, the creeping wave arrivals that travelled around the top cylinderto the bottom cylinder were not detected by the microphone because ofthe attenuation that occurs as waves travel around cylindrical surfaces.82APPENDIX A. TRANSLATION FORMULA FOR THE HANKELFUNCTIONThe following derivation shows how to compute the differential of atranslated Hankel function of the second kind. The Hankel function istranslated from one cylindrical coordinate system to another cylindricalcoordinate system and it is differentiated with respect to the variable p of thenew coordinate system. This method assumes that the argument of theHankel function is much larger than unity so that the method of stationaryphase can be employed to determine a numerical solution. To begin theintegral formula for the Hankel function of the second kind is rewritten here:H(2)(kp) = fd exp( j (kp cos + n-(A. 1)F2it. 3it.where the integration path F travels from ( - Joo) to ( + Joo). Totranslate the Hankel function from one cylindrical coordinate system toanother cylindrical coordinate system as shown in figure A. 1 (taken fromChew [3]) a method presented by Watson [10] is used. To change coordinate83Figure A. 1. Translation of the Hankel Functionsystems p” is replaced by p-p’ I and the following change of variablestaken from Figure A.1where a is defined by the followingIp-p’Isina =p>-p<cos(4-Ø’)(A.2)(A.3)yn Aon0” ip90xI p-p’ I cos a = p< sin(-) (A.4)84The variables p> and p< are defmed as the greater and lesser values of p andp’. When these changes are put back into (A. 1) the following is obtained:H(2)(k p-p’ I) = JdO exp j (k I p-p’ I cos(O+Ø”-x) + n(O+Ø”-x) -(A.5)Thus, the desired integral expression used to translate the Hankel function isobtained by factoring out the exp( jn” ) term and moving it over to the righthand side as shown below:Hn(2)(k I p- p’ I )exp(-jnØ”) = JdO exp j (k I p-p’ I cos(O+p”-a) + n(O-a) -(A.6)If it is assumed that the argument of the Hankel function is large, i.e.k p-p’ I>> 1 then the method of stationary phase can be applied to evaluatethe integral given by equation (A.6). To employ the method of stationaryphase the point at which the phase is not changing needs to be determined. Tofind the stationary phase point the phase term of equation (A.6 ) isdifferentiated with respect to 0 and then the expression is equated to zero andsolved for 0 as shown on the following page:850= (k p - p’ I cos(O+4”-x) + n(O-a))sin(O+”-a)= k I p- p’ I),n<kI p-p’I (A.7)p-p’I (A.8)The next step is to differentiate (A.6 ) with respect to p, and then set p = a,where a = the radius of the cylinder. In carrying out the differentiation ofequation (A.6 ) with respect to p the chain rule is used because all of theangles functions of the variable p. To start, the first phase term isdifferentiated with respect to p. This is accomplished by removing theabsolute value sign and rewriting the first exponential term by substitutingequations (A.3 ) and ( A.4 ) into the first phase term of equation ( A.6) asshown below:(p-p’) cos(O+Ø”-x)=(p-p’) cos(O+Ø”)cosx + (p-p’) sin(O+Ø”)sina=p sin(Ø-Ø’)cos(O+Ø”) + ( PS-P cos(Ø-’)) sin(O-i-”)(A.9)86Differentiating (A.9 ) with respect to p gives:[(p-p’) cos(O+Ø”-a)] = sin(Ø-Ø’)cos(9+Ø”) + p sin(Ø-Ø’) [cos(O+Ø”)]- cos(-’)sin(9+”) + (p’-p cos(Ø-’)) [sin(O+Ø”)](A.1O)Carrying through with the differentiation defined in (A.1O ) results in:[cos(O+”)] = - sin(O+”) (0+0”) (A. 11)[sin(0+Ø”)] = cos(0+Ø”) (0+0”) (A.12)and making use of the stationary phase point given by (A.7 ) and (A.8)yields:)),n<kI p-p” (A.13)),n>kI p-p’ (A.14)87a__________(sin1(kIpflp,I))ap(cosh4(kIppII ))=where1)-1-nk ( I p- p’ I )2a( p- p’)ap(A.17)Differentiating the o’. term of ( A. 13 ) and (A. 14) and utilizing equations(A.3)and(A.4)gives:aa a 1 p> - p< cos (0 -—= —(tan- ( . ))ap ap p<sm(Ø-Ø)-p>sin(O-0’)= (p>)2 + (p<)2-2p< p> cos(0- 0’) (A.15)and differentiating the final term of ( A.13 ) and (A.14 ) completes the stepsin differentiating (A.6):1—(A.16)I p-p’ (p>)2+(p<)22p<p>cos0 (A.18)88anda(p-p’) p<-p>cosØ-J(p>)2+(p<)22p<p>cosØNow that the differentiation steps have been completed the approximation tothe differential of the integral described by equation ( A.6) can be computedusing the method of stationary phase. The method of stationary phase statesthat an integral of the form, I= f d F exp( jf) for large can beapproximated by:1 inlt iCf 0______1 1 f(iv)(Oo) 5 (f”(0o) )2Ie 2 eJ (°) jkf’(0o) 1 j2kf’(Oo) f’(Oo) (f’(0o))2(A.20)where f(90) is given by:f(Oo)=p< sin(Ø-Ø’)cos 00- (p>-p< cos(-Ø’)) sin Oo + (A.21)f(9o), f’(Oo), f”(Oo), and f(iv)(O0)are the first, second, third, and fourthderivatives of f(0) and Go is the stationary phase point.89APPENDIX B. MATLAB SIMULATION SOFfWARE% headerl.m $% This function initializes all of the constants and %% most of the bessel and hankel fucntions used in the %% program one cyl grd.mclear;F I ELDfield_tempfieldindex—zeros (tot_freq,numyoints);—zeros (tot_freq,numyoints);=zeros (num_freq,num_points);—l : num_freq;YoXoahPSdCw_maxwmindelta_wdelta_xdelta_tphilphi21.06;1.25;0.106;0.235;0.155;=2* (ps+h);..2*pi*l5000;.2*pi*1000;.2*pi*50;—0.02;—0.000005;3*pj/2;•pi / 2;% Translation angles, distances and% corresponding derivitives used in% the function HPRIME.m$$$aum_cyl—2;phit—zeros (num_cyl);pJt—zeros (numcyl);dppt_dp —zeros (num_cyl);phitt—zeros (num_cyl);dphittdp—zeros (num_cyl);% translation *1 from 1 -> 2phit(1,2)—pi/2;p_pt(1,2)dp,pt_dp(1,2) ——1;phitt(l,2)——pi/2;dphitt_dp(1, 2)—0;num_points-125; % number of speaker positionsnurn_f req—281; % number of frequency samplestot f req 2048; % number of frequencies in total% distance to speaker from cyl. 1% initial speaker location$ radius of each cylinder% height of microphone above ground% distance to mic. location% seperation between adjacent cyl’s% speed of sound% heighest f req. of speaker% lowest f req. of speaker% freq. steps from 1000—15000 Hz% sample spacing for speaker% time between collected samples% mic. angle wrt cyl.% mic. angle wrt image cyl.% translation *2 from 2—> 1phit(2,1) ——pi/2;p_pt(2,1)dp_pt_dp(2,1) —--1;phitt(2,1)—pi/2;dphittdp(2, 1)’0;— (w_min+(jndex-1) *delta_w) Ic;(1 : tot_f req) *dejtat*c;ceil(max(k)*a)+1;—1 /2*tot_freq*delta_w;% AU. the necessary Hankel %% and Bessel functions are %% pre—calculated into a% matrixindex 1—0 : nmax;index2—0 n_max+1;offsetl—n max+1;offset2—o?fsetl+1;Hnps (index, indexl+offsetl)Hnps_im(index, indexl+offsetl)Hna (index, index2+offset2)for nl——n_max:—l,Hnps (index, n1+offst1)Hnps_ixn(index, nl+offsetl)endfor n2——(n_xnax+1):-1,Hna (index,n2+offset2)end% vector of wave numbers% distance travelled at each delta_t% max I of orders sunutied over% Nyquist frequency—besselj (indexl, (k*ps))_*beseely(index1, (k*ps)’);=besselj (indexi, (k* (ps+2*h))’)_j*ssely (indexi, (k* (ps+2*h))’);=besselj (index2, (k*a)’)..._j*ssely(jndex2, (k*a)’);—(—1) (—ni) *Hnps(jndex,offsetl_nl);—(—1) A (—ni) *Hnps_im (index, offsetl—nl);= (—1) A (—n2) *Hna (index, offset2—n2);% The derivitives of the% Hankel and Bessel functions %% are pre—calculatedfor n3——n_max:n_max;Hnap(index,n3+offsetl)endJnap—reai. (Hnap);=k’ /2. * (Hna (index, n3-l+offset2)—Hna (index, n3+1+offset2));% Load in the source wavelet% in the frequency domain to% weight against the Green’s% functionfid— fopen(’inwave2.bin’,’r’);[inwave count] — fread(fid,23000,’short’);INC— fft(inwave,46000);INWAVE— INC(138:7:2098);krn_maxfn90% Load in the reference% wavelet in the frequency% domain to be used for% matched filtering 91fid2 — fopen(’refwave2.flb’ ,‘r’);(refc count]— fread(fid2,16384,’ushort’);REFC— fft(refc);REFWAVE — REFC(82:4:1202);fclose(’all’);end92Simulation of Diffractionby Cylinders of ThreeDifferent Sizeshorizontal trolleyanaaanaaaasrn===nnnII/ \ < speaker9’9’.9’0 < cylinderI I9’ ground I I <—— microphone 9’9’ 9’9’ 9’!sleep 20000for i—l:ceil(num_points/2),x Xo_(i_l)*delta_x;pgl _(xA2+Y0A2)A0.5;,pg2—(x’2+(Yo+d)2)”0.5;p1— (x’2+ (Yo+ps) 2) ‘0.5;p2 .(x’2+(Yo+ps+2*h)%2Y*0.5;phiol —pi/2—atan (x/Yo);phio2—pif2-atan (x/ (Yo+d));phiolm 2*pi_phio2;phio2rn _2*pi_phiol;Hnll (index,indexl+offsetl) —besselj (indexi, (k*pgl)’)_j*bessely(jndexl, (k*pgl)’);Hn12(index,indexl+offsetl) -besselj(indexl, (k*pg2)’)._j*bessely(indexl, (k*pg2)’);for n4——nmax:—l,Hnll(index,n4+offsetl) _(_l)’i_n4)*Hnhl (index, offset1n4;Mn12(index,ri4+offsetl) s.(_l)’(_n4)*Hnl2(index,offset1_n4);endHol —besselj(0, (k*pl)’)_j*bessely(O, (k*pl)’);Ho2 —besselj(0, (k*p2)’)_j*bessely(O, (k*p2)’);for m=l:nuxn_f req,mninax — ceil(k(m)*a)+l;9’ 9’9’ We first must determine the unknown 9’9’ scattering coefficients for each position 9’9’ frequency and order i.e. the An and Anim $9’ using the scattering matrix T and the 9’9’ incident field given by the vector f 9’$index3——nmax : nmax;offset 3—nmax+ 1;psil — Jnap(m,index3+offsetl)...*Hnhl (m,index3+offsetl)....*exp(_j*index3* (phil—phiol))+ Jnap(m,index3+offsetl)...• *Hn12 Cm, index3+offsetl)...*exp (_j*jndex3* (phil—phiolm) ) ;•psi2 — Jnap(m,index3+offsetl)...• *14j Cm, index3+offseti)...• *exp(_j*index3* (phi2—phio2)).+ Jnap(m,index3+offsetl)...• *Hnhl (m, index3+offseti)• *exp (_j*index3* (phi2-phio2m));alpha —[—pail.’; —psi2.’);Li —zeros (2*nlnax+l);L2 _zeros(2*nlnax+i);L3 —zeros (2*nmax+i);L4 zeros(2*runax+i);Ii HPRIME(d , p_pt(2,i), dp_pt_dp(2,i), k(m), phii,...phit(2,l) ,phitt(2,i) ,dphitt_dp(2,l), nrnax);12 —HPRIME(d, p_pt(i,2), dpyt_dp(i,2), k(m), phi2,...phit(l,2),phitt(l,2), dphitt_dp(i,2), nxnax);for n5——nmax:nmax,Li (n5+offset3,n5+offset3)— Hnap(m,n5+offsetl) *exp(_j*n5*phji);L2 (nS+offset3,n5+offset3)— Ii (n5+nmax+1);L3 (n5+offset3,n5+offset3)— 12 (n5+nmax+i);L4 (n5+offset3,n5+offset3)— Hnap(m,n5+offseti)*exp(_j*n5*phi2);endT—[ Li L2 ;...L3 L4 );f — inv(T)*alpha;Anl _f(l:2*nxnax+i); % scattering coefficient for cyl. iAn2 •f(2*nlnax+2:4*nnlax+2); % scattering coefficient for cyl. 2% Now that the scattering coefficients have %% been determined the field can be% calculated in the frequency domainterml(i:2*ninax+i)Anl.’ .*Hnps(m,index3+offseti)•*exp(_j*jndex3*phjl);term2(l:2*nxnax+i)_An2.’ .*Hnps im(m,index3+offseti)....*exp(_j*index3*phi2);fielci(m, i) —INWAVE (m) *zorjj (REFWAVE Cm) ) *exp(_j*pj/4) *k Cm) ‘O.5*...(Hoi Cm) +1402 Cm) +sum(terml+tertn2));clear psil psi2;clear alpha;clear Li L2 L3 L4;clear Ii 12;clear T;clear f;clear Ani An2;clear termi term2;endclear Hnli Rn12 Hoi Ho2 ;$Before taking the ifft we must $fill in the negative frequencies. $We do this by noting that for $real data the negative frequency %values are the complex conjugate $of the positive frequency values $$for s—O nurn_freq-i,field texnp(((—w max+2*fn)/delta w+i)+s,i)—con(field(num_freq—s,i));field_temp((w_mtn/delta_w+i)+s,I)=field(i+s,i); 94endFIELD(i:tot_freq,i)_(i./(32*pi*r)).’O.5.I.*ifft(field_ten1p(i:tot_freq,j));end% The other half of the speaker% poaitions are filled in by noting %% that the model is sytnetric about the %% x—O plane. $for i2—ceil(nuinj,oints/2)+i:nuxnpoints,FIELD(l:tot_freq,i2)—FIELD(l:tot_freq,numyoints+i—(i2-i));endFIELD1_lg—FIELD;save(’FIELDi_lg’);end% Function HPRIME(pt,p_pt,dp_ptdp,k,phi,phit,phitt,dphitt_dp,nmax) %% This function computes the derivitive of the% hankej. function Hn(kIp_pI)*exp(_*n* (phiphi))% by taking the devitive of the solution obtained% from the integral expression by the method of% steepest descent% The derivitive is taken with respect to p, andI then p is made equal to a (p—a).II Inputs :— distance vectors ( Pt and p_ptI- derivitive of p_pt w.r.t p (dp_pt_dp)I— angles of each distance vectorI C phi, phit, phittI— derivitive of phitt w.r.t. p ( dphitt_dpI— order of Hankel derivitive C nI— wave number beta (Ic)II Outputs:— derivjtjve of Mankel function aboveI w.r.t. p and p is equal to aIfunction (Hnd)—HPRIME(pt,p_pt,dp_pt_dp,k,phi,phit,phitt,dphitt_dp,nlnax)ps—O.O57; I radius of cylinderpgpt; I distance between cylindersif sin (phi—phit)——0,aipha—pi /2;elseif sin(phi—phit)>0alpha—atan( (pg_ps*cos (phi—phit))/ (ps*sin(phi_phit)));else,alpha—pi+atan( (pg_ps*cos(phi_phit)) / (ps*sin(phi_phit)));endlambdaj*k;dalpha_dp=_pg*sin (phi-phit) / (psA2+pgA2_2*ps*pg*cos (phi-phit));for n——nmax: nxnax,if n<k*p_pt,theta— pi + alpha- phitt - asin(n/(k*p_pt));term2 — l/(l_(n/(k*p_pt))’2)0.5 * (_n)/(k*p_pt2) * dp_pt_dp;elsetheta— pi/2 + alpha— phitt - j*acosh(n/(k*p_pt));term2 — l/((n/(k*p_pt))A2_l)A0.5 * (_j*rz)/(k*p_pt’2) * dp_pt_dp;endF — l/pi*exp(_*n*pi/2)*...(j*k* (sin (phi—phit) (theta+phitt)_ps*sin (phi—phit) *j (theta+phitt) * (dalpha_dp—term2)—cos (phi—phit)*3j(theta+phitt).+ (pg_ps*cos (phi-phit) ) *c (theta+phitt) * (dalpha_dp—term2))+j*n* (—dphitt_dp—term2));Fl— l/pi*exp(_j*n*pi/2)**k...* (—sin (phi—phit) *5jfl (theta+phitt)ps* sin (phi—phit) (theta+phitt) * (dalpha_dp—term2)—cos (phi —phit) *cos (theta+phitt)— (pg_ps*cos (phi—phit)) *ji (theta+phitt) * (dalpha_dp—terxn2));F2 l/pi*exp(_j*n*pi/2)**k...* (—sin (phi—phit) *co (theta+phitt)+ps*sin (phi—phit) *gjfl (theta+phitt) * (dalpha_dp—term2)+cos (phi-phit) *3j(theta+phitt)...— (*cos (phi—phit) ) *c (theta+phitt) * (dalpha_dp—terni2));f ps*sjn(phj_phjt)*cos(theta+phjtt)...+ (pg—ps *cos (phi—phit) ) *3jfl (theta+phitt)+ i* (theta—alpha) /k;fl__ps*sin (phi—phit) *sin(theta+phitt)+ (pq_ps*cos (phi—phit) ) *Q5 (theta+phitt)* sin (phi—phit) *cc (theta+phitt)...— (pg_ps*coz (phi—phit) ) *j (theta+phitt);f 3— ps*sin(phi_phit)*sin(theta+phitt)...(pg_ps*c (phi—phit) ) (theta+phitt);f4.—f2;Hnd(n+nmax+l) — (_2*pi/(larnbda*f2)[AO.5 * F * exp(larnbda*f)...* (1+1/ (j*2*k*f2) * C (f3/f2) * (Fl/F)+1/4* (f4/f2) _5/12* (f3/f2) A2_ (F2/F)));endend97REFERENCES[1] V. Twersky, “Multiple Scattering of Radiation by an ArbitraryConfiguration of Parallel Cylinders”, Journal of the Acoustical Society ofAmerica, Vol. 24, P. 42, 1952.[2] V. Twersky, Journal of the Acoustical Society ofAmerica, Vol. 23,p. 407, 1952.[3] Weng Cho Chew, “Waves and Fields in Inhomogeneous Media”,New York: Van Nostrand Reinhold, 1990. pp. 465-469.[4] 3. E. Burke, D. Censor, and V. Twersky, Journal of the AcousticalSociety ofAmerica, Vol. 37, p. 5, 1965.[5] V. E. Glazanov, “Diffraction of a Wave Radiated by a Cylinder at aGrating of Acoustically Compliant Cylinders”, in Soviet Acoustics-Physics, Vol. 14, No. 4, April-June, 1969, pp 447-450, AmericanInstitute of Physics, 1969.[6] R. P. Radlinski and T. 3. Meyers, Journal of the Acoustical SocietyofAmerica, Vol. 56, p. 842, 1974.[7] J. W. Young and 3. C. Bertrand, Journal of the Acoustical SocietyofAmerica, Vol. 58, pp. 1190-1195, 1975.[8] G. 0. Olaofe, Radio Science, Vol. 5, p. 1351, 1970.[9] S. N. Yoursi and F. 3. Fahy, Acustica, Vol. 29, p. 278, 1973.[10] G. N. Watson, “A Treatise on the Theory ofBessel Functions”,pp. 358 - 359, Cambridge: University Press, 1945.98[11] Jin Au Kong, “Electromagnetic Wave Theory”, second edition,pp. 304 - 311, Toronto, John Wiley & Sons Inc., 1990.[12] Garfield Richard Mellema, “An Acoustic Scatter-Mapping ImagingSystem”, M. A. Sc. Thesis. U. B. C., Department of ElectricalEngineering, 1990.[13] Constantine A. Balanis, “Advanced Engineering Electromagnetics”,pp. 626-633, John Wiley & Sons, 1989.[14] M. A. Abramowitz and I. A. Stegun, eds., “Handbook ofMathematical Functions”, Ch. 9, Sec. 9.1.79, p. 363, Washington, D.C.: National Bureau of Standards, 1964.[15] J. J. Bowman, T. B. A. Senior, and P. L. E. Uslenghi (Eds.),“Electromagnetic and Acoustic Scattering by Simple Shapes”, pp. 115 -116, Amsterdam: North-Holland, 1969.


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items