Source Code for Module wiat.MM.Waveguides

  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  """ Waveguides.py - Definition of waveguides for use with the multimodal 
 23      wave propagation solver. 
 24   
 25  """ 
 26   
 27  import math 
 28  import numpy 
 29  from numpy import sin, cos, cosh, arcsin, sqrt, real, imag, pi, log10 
 30   
 31 -class WaveguideProfile: 
 32 -    def __init__(self): 
 33          print "test" 
 34   
 35 -    def assertZ(self,z): 
 36          pass 
 37          #assert 0. <= z <= self.L, "position out of range" 
 38   
 39 -    def get_length(self): 
 40          return self.L 
 41   
 42 -    def radius(self,z): 
 43          " to be defined be derived class " 
 44          pass 
 45   
 46 -    def radiusp(self,z): 
 47          " to be defined by derived class " 
 48   
 49 -    def S(self,z): 
 50          " return the area of the axisymmetric duct from free end to source " 
 51          self.assertZ(z) 
 52          return numpy.pi * self.radius(z)**2 
 53   
 54 -    def Sp(self,z): 
 55          self.assertZ(z) 
 56          return 2 * numpy.pi * self.radius(z) * self.radiusp(z) 
 57   
 58 -    def Sp_S(self,z): 
 59          self.assertZ(z) 
 60          return 2 * self.radiusp(z) / self.radius(z) 
 61       
 62 -class Cone(WaveguideProfile): 
 63 -    def __init__(self,d0,de,L): 
 64          self.a0 = d0/2 
 65          self.ae = de/2 
 66          self.L = L 
 67 -    def radius(self,z): 
 68          self.assertZ(z) 
 69          return self.ae + (self.a0-self.ae)*z/self.L 
 70       
 71 -    def radiusp(self,z): 
 72          self.assertZ(z) 
 73          return (self.a0-self.ae)/self.L 
 74   
 75 -class Cylinder(WaveguideProfile): 
 76 -    def __init__(self,d,L): 
 77          self.a = d/2 
 78          self.L = L 
 79 -    def radius(self,z): 
 80          self.assertZ(z) 
 81          return self.a 
 82       
 83 -    def radiusp(self,z): 
 84          self.assertZ(z) 
 85          return 0. 
 86   
 87 -class QuarterCircleHorn(WaveguideProfile): 
 88 -    def __init__(self,d0,de,theta1=0.,theta2=numpy.pi/2.): 
 89          """ 
 90          d1: input diameter of the horn 
 91          d2: output diameter of the horn 
 92          theta1: input angle of the horn (radians) 
 93          theta2: output angle of the horn (radians) 
 94   
 95          a: input radius of the horn 
 96          r: radius of the quarter circle 
 97          """ 
 98          self.d0 = d0 
 99          self.de = de 
100          self.theta1 = theta1 
101          self.theta2 = theta2 
102   
103          self.a0 = d0/2. 
104          self.r_horn = (de-d0) / (2. * (cos(theta1)-cos(theta2))) 
105          self.L = self.r_horn * (sin(theta2)-sin(theta1)) 
106  #        if theta2 == numpy.pi/2.: 
107  #            theta2 = theta2 - 0.000001  
108   
109          self.ze = self.r_horn * sin(theta2) 
110          self.y0 = self.r_horn * cos(theta1)  
111   
112 -    def __str__(self): 
113          return "QuarterCircleHorn(d0=%.4f,de=%.4f,theta1=%.2f,theta2=%.2f)"%\ 
114                  (self.d0,self.de,self.theta1,self.theta2) 
115 -    def radius(self,z): 
116          self.assertZ(z) 
117          theta = arcsin( (self.ze - z) / self.r_horn) 
118          return self.a0 + self.y0 - self.r_horn * cos(theta) 
119   
120 -    def radiusp(self,z): 
121          """ derivative of the radius along z 
122              the 0.000001 value was added to avoid infinite derivatives 
123          """ 
124          self.assertZ(z) 
125          arg = (self.ze - z) / self.r_horn 
126          theta = arcsin( arg ) 
127          if (1. - arg) < 1e-15: 
128              arg = 1. - 1.e-15 
129          dtheta = -1. / (self.r_horn * sqrt(1.-arg**2)) 
130          return  self.r_horn * sin(theta) * dtheta 
131   
132 -class BesselHorn(WaveguideProfile): 
133 -    def __init__(self,d1,d2,L,alpha): 
134          self.L = L 
135          self.alpha = alpha 
136          self.B = 10**(log10(L) - log10(d1**-1-d2**-1)) 
137          self.z0 = 10**( alpha**-1 * log10(self.B/d2) ) 
138   
139 -    def radius(self,z): 
140          return 0.5 * self.B / (self.L - z + self.z0)**self.alpha 
141   
142 -    def radiusp(self,z): 
143          return 0.5 * self.B * self.alpha / \ 
144                  (self.L - z + self.z0)**(self.alpha+1.) 
145   
146  #class ExponentialHorn(WaveguideProfile): 
147  #    def __init__(self,d1,d2,L,theta1): 
148  #        self.L = L 
149  #        self.a1 = d1/2. 
150  #       self.a2 = d2/2. 
151  #        self.m = numpy.tan(math.radians(theta1/2.)) / self.a1 
152  # 
153  #    def radius(self,z): 
154  #        return numpy.sqrt( self.S0 * numpy.exp(self.m * z) / numpy.pi ) 
155  # 
156  #    def radiusp(self,z): 
157  #        pass 
158  # 
159