Source Code for Module wiat.MM.Radiation

  1  # Copyright 2008 Antoine Lefebvre 
  2  # 
  3  # This file is part of "The Wind Instrument Acoustic Toolkit". 
  4  # 
  5  # "The Wind Instrument Acoustic Toolkit" is free software: 
  6  # you can redistribute it and/or modify it under the terms of 
  7  # the GNU General Public License as published by the 
  8  # Free Software Foundation, either version 3 of the License, or 
  9  # (at your option) any later version. 
 10  # 
 11  # "The Wind Instrument Acoustic Toolkit" is distributed in the hope 
 12  # that it will be useful, but WITHOUT ANY WARRANTY; without even 
 13  # the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
 14  # See the GNU General Public License for more details. 
 15  # 
 16  # You should have received a copy of the GNU General Public License 
 17  # along with "The Wind Instrument Acoustic Toolkit". 
 18  # If not, see <http://www.gnu.org/licenses/>. 
 19  # 
 20  # All rights reserved. 
 21   
 22  """ Radiation.py - Multimodal radiation matrix implementations. 
 23   
 24  """ 
 25   
 26  import os 
 27  import numpy 
 28  from numpy import sin, cosh, sqrt, real, imag, pi 
 29  import scipy 
 30  import scipy.integrate 
 31  import scipy.interpolate 
 32  from scipy.interpolate import UnivariateSpline 
 33  import scipy.io 
 34  import scipy.special 
 35  from scipy.special import j0,j1,jn_zeros,jvp 
 36   
 37  from Functions import flat2matrix, matrix2flat 
 38   
 39 -class RadiationFlanged: 
 40      """ Class implementing the Zorumski multimodal radiation impedance matrix. 
 41       
 42          [1] N. Amir, V. Pagneux, and J. Kergomard, 
 43          "A study of wave propagation in varying cross-section waveguides by  
 44          modal decomposition. Part II. Results," 
 45          The Journal of the Acoustical Society of America,   
 46          vol. 101, May. 1997, pp. 2504-2517. 
 47   
 48          [2] W.E. Zorumski, "Generalized radiation impedances and reflection  
 49          coefficients of circular and annular ducts," 
 50          The Journal of the Acoustical Society of America,   
 51          vol. 54, Dec. 1973, pp. 1667-1673. 
 52   
 53          [3] Jonathan Kemp, "Theoretical and experimental study of wave  
 54          propagation in brass musical instruments", 
 55          Ph.D. Thesis, University of Edinburgh, 2002         
 56           
 57      """ 
 58      N = 15 
 59       
 60 -    def __init__(self,diameter,number_of_modes = N, do_interpolation = True): 
 61          self.radius = diameter / 2. 
 62          self.number_of_modes = number_of_modes 
 63   
 64          assert number_of_modes >=1, 'Minimal number of modes is 1' 
 65          if number_of_modes == 1: 
 66              self.gamma = [0.] 
 67          else: 
 68              self.gamma = numpy.concatenate( 
 69                      ([0.],jn_zeros(1,self.number_of_modes-1))) 
 70   
 71          # if there is a data file, load it and do interpolation 
 72          self.do_interpolation = do_interpolation 
 73          self.backup_filename = os.path.join(os.path.dirname(__file__),'Znm') 
 74          if do_interpolation: 
 75              try: 
 76                  self.data_table = scipy.io.loadmat(self.backup_filename) 
 77                  N = int(self.data_table['number_of_modes']) 
 78                  Znm = numpy.reshape(self.data_table['Z_nm'], 
 79                                      (N,N,self.data_table['length_ka'])) 
 80                  self.data_table['Z_nm'] = Znm 
 81                  #print 'Creating interpolation splines for the multimodal  
 82                  #radiation.' 
 83                  self.spl_real = [] 
 84                  self.spl_imag = [] 
 85                  for i in range(N): 
 86                      self.spl_real.append([]) 
 87                      self.spl_imag.append([]) 
 88                      for j in range(N): 
 89                          #print self.data_table['ka'] 
 90                          #print numpy.real(self.data_table['Z_nm'][i][j]) 
 91                          self.spl_real[i].append( 
 92                                  UnivariateSpline(self.data_table['ka'], 
 93                                      numpy.real(self.data_table['Z_nm'][i][j]), 
 94                                      s=0.0000001)) 
 95                          self.spl_imag[i].append( 
 96                                  UnivariateSpline(self.data_table['ka'], 
 97                                      numpy.imag(self.data_table['Z_nm'][i][j]), 
 98                                      s=0.0000001)) 
 99               
100              except IOError: 
101                  #print arg 
102                  print 'Calculating radiation data table for interpolation.' 
103                  # if there is no data file, create it 
104                  ka = numpy.arange(0.,15.,0.1) 
105                  self.do_interpolation = False 
106                  self.data_table = {'number_of_modes': self.number_of_modes, 
107                                     'length_ka': numpy.size(ka,-1), 
108                                     'ka': ka, 
109                                     'Z_nm': self.calculateZnm(ka/self.radius)} 
110                  self.do_interpolation = True 
111                  scipy.io.savemat(self.backup_filename, self.data_table) 
112           
113 -    def __call__(self,N,k): 
114          assert N <= self.number_of_modes,\ 
115              'number of modes exceed radiation impedance maximum' 
116          Z = numpy.empty((N,N),dtype=complex) 
117          for i in range(N): 
118              for j in range(N): 
119                      Z[i][j] = self.Z_ij(i,j,k) 
120          return matrix2flat(Z) 
121       
122 -    def Phi(self,r,i): 
123          return j0(self.gamma[i]*r/self.radius)/j0(self.gamma[i]) 
124   
125 -    def D(self,tau,i,k): 
126          """ from [3] eq. 4.1 """ 
127          lambdan = self.gamma[i]/(k*self.radius) 
128          if i == 0 and tau == 0: return 0. 
129          return -sqrt(2.)*tau*j1(tau*k*self.radius) / (lambdan**2-tau**2) 
130   
131 -    def integrand_Z_ij_real(self,phi,i,j,k): 
132          return sin(phi) * self.D(sin(phi),i,k) * self.D(sin(phi),j,k) 
133   
134 -    def integrand_Z_ij_imag(self,xhi,i,j,k): 
135          return cosh(xhi) * self.D(cosh(xhi),i,k) * self.D(cosh(xhi),j,k) 
136               
137 -    def Z_ij(self,i,j,k): 
138          #print i,j 
139          if self.do_interpolation: 
140              ka = k * self.radius  
141              if ka >= self.data_table['ka'][-1]: 
142                  print 'Warning: value of ka not in radiation matrix table' 
143              Z_ij = complex(self.spl_real[i][j](ka)[0], 
144                             self.spl_imag[i][j](ka)[0]) 
145          else: 
146              if k == 0.: 
147                  Z_ij = complex(0.,0.) 
148              else: 
149                  realpart = scipy.integrate.quad(self.integrand_Z_ij_real, 
150                          0., numpy.pi/2,args=(i,j,k),limit=5000) 
151                  imagpart = scipy.integrate.quad(self.integrand_Z_ij_imag,  
152                          0., 100, args=(i,j,k),limit=5000) 
153   
154                  # TODO- Check the sign of the imaginary part  
155                  # (kemp is positive, Zorumski negative) 
156                  # Antoine Lefebvre (2008) - It seems that the negative sign  
157                  # is correct, 
158                  # from comparison with plane wave models 
159                  Z_ij =  complex(realpart[0], -imagpart[0]) 
160   
161          return Z_ij 
162       
163 -    def calculateZnm(self,ks,number_of_modes=N):         
164          Zs = numpy.empty((number_of_modes, number_of_modes,len(ks)), 
165                           dtype=complex) 
166          l = 0 
167          print "Calculating flanged modal radiation matrix..." 
168          for k in ks: 
169              for i in range(number_of_modes): 
170                  for j in range(number_of_modes): 
171                      Zs[i][j][l] = self.Z_ij(i,j,k) 
172              l = l+1 
173          print "done." 
174   
175          return Zs 
176           
177   
178   
179  # this is not correct 
180  # a formulation for the multimodal radiation impedance have to be found 
181  # maybe in 
182  # 
183  # Wang, K. S., & Tszeng, T. C. (1984). 
184  # Propagation and radiation of sound in a finite length duct. 
185  # Journal of Sound and Vibration, 93(1), 57-79. 
186 -class RadiationUnflanged: 
187 -    def __init__(self,radius): 
188          self.a = radius 
189   
190 -    def __call__(self,N,k): 
191          ka = k * self.a 
192          R0 = (1 + 0.2*ka - 0.084*ka**2)/(1 + 0.2*ka + (0.5 - 0.084)*ka**2) 
193          delta = 0.6113 * ((1+0.044*ka**2)/(1+0.19*ka**2) - 0.02*(numpy.sin(2*ka))**2 ) 
194          R = - R0 * numpy.exp( -2 * 1j * ka * delta) 
195          Zl = (1+R)/(1-R) 
196          M = numpy.diag(Zl*numpy.ones(N)) 
197          return matrix2flat(M) 
198   
199       
200   
201 -class RadiationIdealOpen: 
202      pass 
203   
204 -class CharacteristicImpedance: 
205      pass 
206   
207 -def AdmittanceClosed(N,k): 
208      """ Rigid wall admittance matrix is zero (N x N) 
209          they are complex valued so we should return a flat vector of 2*N**2 values 
210      """ 
211      return numpy.zeros(2*N**2) 
212   
213  #    Z = numpy.diag( self.Zc(0.,k)/10. , 0) 
214  #        return numpy.concatenate( (numpy.real(Z.flatten()), numpy.imag(Z.flatten()) )) 
215  #    pass 
216