##----------------------------------------------------------------------------
## Name: RotarySpeaker.py
##
## Compute the angular pressure (beampattern) of the rotary speaker.
## The wall vibrations can be included with -v & -W
##
## Coordinates are: x - horizontal (x=0 along wall/baffle)
## y - vertical (y=0 along wall/baffle)
## z - fan axis
## The beampattern is taken in the x-z plane
## Origin (0,0) is the fan center
##
## Depends on Python SciPy and NumPy
##
## Author: J Park
##
## Created: 4/9/08
##
## Modified:
##----------------------------------------------------------------------------
import os, sys, string, getopt, time
Version = sys.argv[0]+" : "+str( time.asctime( ( time.gmtime(time.time()) ) ) )
DEBUG = False
# Python array module and numpy/scipy are NOT compatible as standalone modules.
# numpy/scipy requires that the array constructor must be passed an object
# that is either an array, or list: a = array([]) will create an empty array.
#
# With the python array module, one can define an array object as:
# a = array('d'), and then call a.append(1.2) to build the array as
# needed. In numpy/scipy, you have to first create an array object a,
# then call b = append(a, 1.2) to return an array b that appends the
# value 1.2 to the existing array object a.
#
try:
from numpy import * # www.numpy.org numpy.scipy.org
except ImportError:
print "Failed to import numpy."
try:
import pylab as mp # matplotlib.sourceforge.net
except ImportError:
print "Failed to import pylab."
try:
from scipy import special as sp # www.scipy.org
except ImportError:
print "Failed to import scipy."
# Global variables
# Default Parameters
PlotNum = 0 # -p
WriteOutput = False # -w
Monopole = False # -m compute monopole field only
FlexBaffle = False # -v include baffle
AxisOnly = False # -z
OutputFile = "piston.dat" # -o
Frequency = 16. # -f blade actuation freq (Hz)
R_i = 4. * 0.0254 # -a piston inner radius (in) to (m)
R_o = 10. * 0.0254 # -b piston outer radius (in) to (m)
dL = 1.0 * 0.0254 # -l piston grid resolution (in) to (m)
F_r = 10. # -R Field range (m)
range_inc = 0. # -r recurring range calculation increment (m)
B_theta_max = 35. * pi/180. # -x blade max actuation angle (deg) to (rad)
# this is set by Amp Gain and Input Voltage (RMS)
#Frequency_fan= 11.02 # -F fan rotational speed (Hz) Cfan=23, 661 RPM
Frequency_fan= 16.8 # -F fan rotational speed (Hz) Cfan=35, 1006 RPM
delta_wall = 0.003 # -W wall displacement (m)
c_0 = 343. # sound speed (m/s)
rho = 1.2 # air density (kg/m^3)
Omega = 2.*pi*Frequency # (rad/s)
Omega_fan = 2.*pi*Frequency_fan # fan rotational speed (rad/s)
K = Omega / c_0 # wavenumber
T = 0. # time
A_total = pi*(pow(R_o,2) - pow(R_i,2)) # area of 'piston' (m^2)
V_r = Omega_fan * (R_i + (2.*(R_o - R_i)/3.)) # blade eff forward velocity (m/s)
V_r_y = V_r * sin( B_theta_max ) # forward velocity normal to blade
B_w = 2. * 0.0254 # blade width from actuator axis (m)
V_b = Omega * 2.* B_w / 3.# blade eff actuation velocity normal to blade (m/s)
V_b_avg = V_r_y + V_b # effective particle velocity normal to blade (m/s)
C_fan = ( rho * c_0 * V_b_avg * K ) / ( 2. * pi )
V_wall = 2. * delta_wall * Frequency # wall velocity (m/s)
C_wall = ( rho * c_0 * V_wall * K ) / ( 2. * pi )
F_theta = arange( 0., pi/2., pi/180. ) # Field angles (rad)
wall_dL = 0.25
wallHeightMin = -0.3
wallHeightMax = 3.2
wallWidthMin = -5
wallWidthMax = 3.8
wallHeight = wallHeightMax - wallHeightMin
wallWidth = wallWidthMax - wallWidthMin
bldg_f0 = c_0/30. # c/lambda_0 = 11.43 Hz, lambda_0 = 2 L
bldg_omega0 = 2. * pi * bldg_f0
bldg_k0 = bldg_omega0 / c_0
# List of field points, list preserves order
Field = []
# array of field point pressures for output
P_n_amp = array([], dtype=float32)
P_n_real = array([], dtype=float32)
P_n_imag = array([], dtype=float32)
fanAxisPres = complex(0., 0.)
# Dictionary of piston grid points
# key is piston point (x,y) tuple, value is class object
PistonGrid = {}
WallGrid = {}
# Dictionary of rays from piston grid points to field points
# key is tuple of source/field (PistonPoint, FieldPoint), value is class object
Rays = {}
##----------------------------------------------------------------------------
class PistonGridPoint:
def __init__( self, x, y ):
self.x = x
self.y = y
self.z = 0.
self.A = 0.
self.C = C_fan
PistonGrid[(x,y)] = self
##----------------------------------------------------------------------------
class WallGridPoint:
def __init__( self, x, y ):
self.x = x
self.y = y
self.z = 0.
self.A = 0.
#self.C = C_wall
# sin distribution of velocities, nodes at bottom/top
self.C = C_wall * sin( pi * (self.y - wallHeightMin) / wallHeight ) * \
sin( pi * (self.x - wallWidthMin) / wallWidth)
WallGrid[(x,y)] = self
##----------------------------------------------------------------------------
class FieldPoint:
def __init__( self, r, theta ):
self.r = r
self.theta = theta
#self.y = 0.
self.y = -0.34 # the fan and microphone are not planar
self.x = r * sin( theta )
self.z = r * cos( theta )
self.pressure = complex(0., 0.)
self.rays = []
# insert object into Field container
Field.append(self)
def SumPressure( self ):
global fanAxisPres
# sum the amplitude of this field point from all rays
# this is complex addition of ray amp/phase
# exp[j(kr+wt)] = exp[j phi] = cos(phi) + j sin(phi)
for ray in self.rays:
# Scale for element area, piston/wall velocity and range decay
C = ( ray.gridPoint.A * ray.gridPoint.C / self.r )
self.pressure = self.pressure + C * ( \
complex( cos(ray.phi), sin(ray.phi) ) )
## # This adds a plane wave from the wall
## self.pressure = self.pressure + 0.5 * C * ( \
## complex( cos(ray.phi), sin(ray.phi) ) +
## complex( cos(K*self.z), sin(K*self.z) ) )
if FlexBaffle:
# compute fan pressure alone on the axis
if self.theta == 0.:
if isinstance(ray.gridPoint, PistonGridPoint):
fanAxisPres = fanAxisPres + C * ( \
complex( cos(ray.phi), sin(ray.phi) ) )
if self.theta == 0.:
if not FlexBaffle:
fanAxisPres = self.pressure
print "fan axis pressure", 20.*log10(abs(fanAxisPres))
##----------------------------------------------------------------------------
class Ray:
def __init__( self, gp, fp ):
self.gridPoint = gp
self.fieldPoint = fp
self.r = self.distance(gp, fp)
self.amp = 1./(pow(self.r,2))
self.loss = 20 * log10(self.r)
if isinstance(self.gridPoint, PistonGridPoint):
# propagation from the fan
self.phi = Omega * T - K * self.r
else:
# wall points are out of phase with piston
self.phi = bldg_omega0 * T - K * self.r - bldg_k0 * self.r # + \
#sqrt(pow(K * gp.x, 2) + pow(K * gp.y, 2)) # + 180 * pi/180.
# save reference to this ray in the field point
fp.rays.append( self )
# add ray to main Rays container
Rays[(gp,fp)] = self
def distance( self, gp, fp ):
dx = gp.x - fp.x
dy = gp.y - fp.y
dz = gp.z - fp.z
r = sqrt( pow(dx,2) + pow(dy,2) + pow(dz,2) )
return r
##----------------------------------------------------------------------------
## Main module
def main():
print Version
global fdOutFile, P_n_amp, P_n_real, P_n_imag, F_r, Field, Rays, fanAxisPres
ParseCmdLine()
if DEBUG:
print "K R_o=%.2f" % ( K*R_o ),
print "V_r=%.2f(m/s) V_r_y=%.2f(m/s) V_b=%.2f(m/s) V_b_avg=%.2f(m/s)" % \
(V_r, V_r_y, V_b, V_b_avg)
print "A_total=%.3f" % (A_total),
print "C_fan=%.2f" % (C_fan),
print "C_wall=%.2f" % (C_wall)
# Create a grid of points for the piston
if Monopole:
PistonGridPoint( 0., 0. )
else:
for x in arange( -2*R_o, 2*R_o, dL ):
for y in arange( -2*R_o, 2*R_o, dL ):
xsqr = pow(x,2)
ysqr = pow(y,2)
if xsqr == 0. and ysqr == 0.:
continue
if (xsqr + ysqr < R_o) and (xsqr + ysqr > R_i):
PistonGridPoint( x, y )
# area of each grid point
numPistonGridPoints = len(PistonGrid)
for pg in PistonGrid.values():
pg.A = A_total / numPistonGridPoints
print "%d piston grid points" % (numPistonGridPoints),
# Create a grid for the wall
numWallGridPoints = 0
if FlexBaffle:
for x in arange( wallWidthMin, wallWidthMax, wall_dL ):
for y in arange( wallHeightMin, wallHeightMax, wall_dL ):
xsqr = pow(x,2)
ysqr = pow(y,2)
if (xsqr + ysqr < R_o) :
continue
WallGridPoint( x, y )
# area of each grid point
numWallGridPoints = len(WallGrid)
for wg in WallGrid.values():
wg.A = wall_dL * wall_dL
print "%d wall grid points" % (numWallGridPoints)
if range_inc:
stop = 3
else:
stop = 1
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
for range_num in range(stop):
# reset the containers if multiple ranges
if range_inc:
Field = []
Rays = {}
P_n_amp = array([], dtype=float32)
P_n_real = array([], dtype=float32)
P_n_imag = array([], dtype=float32)
fanAxisPres = complex(0., 0.)
# Create the field points
if DEBUG:
print 'Range', F_r,
for angle in F_theta:
FieldPoint( F_r, angle )
if AxisOnly:
if angle == 0:
break
numRays = len(Field) * (numPistonGridPoints + numWallGridPoints)
print "%d field points," % (len(Field)),
if numRays > 1.2E6:
print "Too many rays:", numRays
sys.exit(1)
# From each piston grid point, compute rays to the field points
for pg in PistonGrid.values():
for fp in Field:
Ray( pg, fp )
print "%d piston rays" % (len(Rays)),
# From each wall grid point, compute rays to the field points
for wg in WallGrid.values():
for fp in Field:
Ray( wg, fp )
print "%d total rays." % (len(Rays))
# At each field point sum the ray pressures
for fp in Field:
fp.SumPressure()
# Extract field results into arrays for plotting
for fp in Field:
#P_n_amp = append(P_n_amp, abs(fp.pressure))
# normalize to the fan axis pressure
P_n_amp = append(P_n_amp, abs(fp.pressure) *
(abs(fanAxisPres) / abs(Field[0].pressure)))
P_n_real = append(P_n_real, fp.pressure.real)
P_n_imag = append(P_n_imag, fp.pressure.imag)
## # Normalize amplitude
## P_min = min(P_n_amp)
## P_max = max(P_n_amp)
## if P_max > P_min:
## P_n_amp = (P_n_amp - P_min)/(P_max - P_min)
## for fp in Field:
## fp.amp = (fp.amp - P_min)/(P_max - P_min)
## # compute ideal directivity for piston in rigid infinite baffle
## P_J1 = array([], dtype=float32)
## for fp in Field:
## arg_j = K * R_o * sin( fp.theta )
## if fp.theta != 0.:
## j1 = 2. * sp.j1( arg_j ) / arg_j
## else:
## j1 = 1.
## P_J1 = append( P_J1, j1 )
## mp.figure(1)
## mp.subplot(2,1,2)
## mp.plot ( F_theta*(180./pi), P_J1, linewidth=2 )
#==========================================================================
# Plotting & Output
#
freqLabel = "Frequency %.0f Hz" % (Frequency)
R_oLabel = "R_o=%.2f in" % (R_o / 0.0254)
R_iLabel = "R_i=%.2f in" % (R_i / 0.0254)
rangeLabel = "Range=%.1f m" % (F_r) # for legend() w/ multiple ranges
#=========================================================
if PlotNum > 5: # plot the field points
x = array([], dtype=float32)
z = array([], dtype=float32)
for fp in Field:
#print fp.x, fp.z
x = append(x, fp.x)
z = append(z, fp.z)
mp.figure(6)
mp.plot ( x, z, '.' )
#=========================================================
# On-axis pressure for a piston in baffle Kinsler & Frey pg 176
if PlotNum > 4:
r_axis = array([], dtype=float32)
for r in arange( 1E-3, 2 * F_r, 2 * F_r / 100. ):
r_axis = append( r_axis, r )
P_axis = abs( sin( (K*r_axis/2) * sqrt( 1 + pow((R_o/r_axis),2)) - 1 ) )
mp.figure(5)
mp.title("On-axis pressure for a piston in baffle")
mp.xlabel("Range (m)", size='16' )
mp.ylabel("Relative Pressure", size='16')
mp.figtext(0.7, 0.25, freqLabel, size='16')
mp.figtext(0.7, 0.2, R_oLabel, size='16')
mp.plot( r_axis, P_axis, linewidth=2 )
#=========================================================
if PlotNum > 3: # real and imaginary pressure vs. angle
mp.figure(4)
mp.subplot(2,1,1)
mp.title("Pressure vs. angle")
mp.ylabel("Real Pressure re Pa", size='16')
mp.figtext(0.6, 0.75, freqLabel, size='16')
mp.figtext(0.6, 0.7, R_oLabel, size='16')
mp.figtext(0.6, 0.65, R_iLabel , size='16')
mp.figtext(0.6, 0.6, rangeLabel, size='16')
mp.plot ( F_theta*(180./pi), P_n_real, linewidth=2 )
mp.subplot(2,1,2)
mp.xlabel("Angle (deg)", size='16' )
mp.ylabel("Imag Pressure re Pa", size='16' )
mp.plot ( F_theta*(180./pi), P_n_imag, linewidth=2 )
#=========================================================
if PlotNum > 2: # plot the piston grid
xp = array([], dtype=float32)
yp = array([], dtype=float32)
xw = array([], dtype=float32)
yw = array([], dtype=float32)
for pg in PistonGrid.values():
#print pg.x, pg.y
xp = append(xp, pg.x)
yp = append(yp, pg.y)
for wg in WallGrid.values():
#print pg.x, pg.y
xw = append(xw, wg.x)
yw = append(yw, wg.y)
mp.figure(3)
mp.plot ( xp, yp, 'r.' )
mp.hold(True)
mp.plot ( xw, yw, 'b.' )
#=========================================================
# polar plot field pressure
if PlotNum > 1:
mp.figure(2)
# This is cool, in numpy.array [::-1] is an 'extended' slice that reverses
polarAngles = concatenate( (-F_theta[::-1], F_theta[1:]) )
polarAmp = concatenate( (P_n_amp[::-1], P_n_amp[1:]) )
mp.polar ( polarAngles, polarAmp, linewidth=2, label=rangeLabel )
#mp.polar ( F_theta, P_n_amp, -F_theta, P_n_amp, linewidth=2, label=rangeLabel )
mp.title("Normalized Beampattern")
mp.figtext(0.28, 0.7, freqLabel, size='16')
mp.figtext(0.28, 0.65, R_oLabel, size='16')
mp.figtext(0.28, 0.6, R_iLabel, size='16')
if not range_inc:
mp.figtext(0.28, 0.55, rangeLabel, size='16')
#=========================================================
# line plot field pressure vs theta
if PlotNum > 0:
mp.figure(1)
#mp.subplot(2,1,1)
mp.title("Pressure Amplitude vs. angle")
mp.xlabel("Angle (deg)", size='16' )
mp.ylabel("Pressure dB re Pa", size='16' )
mp.figtext(0.6, 0.45, freqLabel, size='16' )
mp.figtext(0.6, 0.4, R_oLabel, size='16' )
mp.figtext(0.6, 0.35, R_iLabel, size='16' )
if not range_inc:
mp.figtext(0.6, 0.3, rangeLabel, size='16' )
mp.plot ( F_theta*(180./pi), 20.*log10(P_n_amp), linewidth=2, label=rangeLabel )
# increment the field range for the next run
F_r = F_r + range_inc
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if WriteOutput:
fd = OpenFile( OutputFile, 'w' )
fd.write( "Theta(rad)\tAmp\n" )
for fp in Field:
fd.write( "%.3f\t%.3f\n" % ( fp.theta, fp.pressure ) )
fd.close()
print "Wrote ", len(Field), " points to ", OutputFile
if PlotNum > 0:
# Write out graphic of the plot
plotFileName = OutputFile + ".png"
mp.savefig( plotFileName )
print "Wrote png image to ", plotFileName
if range_inc:
mp.figure(1)
mp.legend( loc='lower left' )
if PlotNum > 1:
mp.figure(2)
mp.legend( loc='lower left' )
# Matplotlib show for all figures (must be only one call)
mp.show()
print "Normal Exit"
## Main module
##----------------------------------------------------------------------------
##----------------------------------------------------------------------------
def usage():
print "Usage: ", sys.argv[0]
print "\t -h (help)"
print "\t -d (debug)"
print "\t -p # plots"
print "\t -w (WriteOutput)"
print "\t -o OutputFile"
print "\t -v FlexBaffle"
print "\t -z AxisOnly"
print "\t -f Frequency (Hz)"
print "\t -F Fan Rotation Frequency (Hz)"
print "\t -a Inner Radius (in)"
print "\t -b Outer Radius (in)"
print "\t -l Grid Resolution (in)"
print "\t -R Field range (m)"
print "\t -r range increment (m)"
print "\t -W wall displacement (m)"
print "\t -x blade actuation angle (deg)"
print "\t -m (monopole)"
##----------------------------------------------------------------------------
def ParseCmdLine():
##
global DEBUG
global delta_wall, C_wall, V_wall, B_theta_max, V_r, V_r_y, V_b, V_b_avg
global Frequency_fan, Omega_fan, C_fan
global PlotNum, WriteOutput, OutputFile, Monopole, FlexBaffle, AxisOnly
global Frequency, R_i, R_o, dL, F_r, Omega, K, range_inc
try:
## Parse the command line arguments in sys.argv via getopt()
opts, args = getopt.getopt(sys.argv[1:], "hdp:wo:f:a:b:l:R:mr:vzW:x:F:")
except getopt.GetoptError, error:
print "ParseCmdLine(): GetoptError: ",error
return
if len(opts) == 0:
usage()
return
argCount = 0
for optElement in opts:
argCount = argCount + 1
if optElement[0] == "-h":
usage()
sys.exit(1)
if optElement[0] == "-d":
DEBUG = True
if DEBUG and 0:
print "opts element %d:" % (argCount), optElement
if optElement[0] == "-p":
PlotNum = int( optElement[1] )
if optElement[0] == "-w":
WriteOutput = True
if optElement[0] == "-o":
OutputFile = optElement[1]
WriteOutput = True
if optElement[0] == "-f":
Frequency = float( optElement[1] )
Omega = 2.*pi*Frequency
K = Omega / c_0
V_wall = 2. * delta_wall * Frequency
V_b = Omega * B_w / 2.
V_b_avg = V_r_y + V_b
C_fan = ( rho * c_0 * V_b_avg * K ) / ( 2. * pi )
C_wall = ( rho * c_0 * V_wall * K ) / ( 2. * pi )
if optElement[0] == "-F":
Frequency_fan = float( optElement[1] )
Omega_fan = 2.*pi*Frequency_fan
V_r = Omega_fan * (R_o + (R_i - R_o/2.))
V_r_y = V_r * sin( B_theta_max )
V_b_avg = V_r_y + V_b
C_fan = ( rho * c_0 * V_b_avg * K ) / ( 2. * pi )
if optElement[0] == "-a":
R_i = float( optElement[1] ) * 0.0254
if optElement[0] == "-b":
R_o = float( optElement[1] ) * 0.0254
if optElement[0] == "-l":
dL = float( optElement[1] ) * 0.0254
if optElement[0] == "-R":
F_r = float( optElement[1] )
if optElement[0] == "-r":
range_inc = float( optElement[1] )
if optElement[0] == "-m":
Monopole = True
if optElement[0] == "-v":
FlexBaffle = True
if optElement[0] == "-W":
delta_wall = float( optElement[1] )
V_wall = 2. * delta_wall * Frequency
C_wall = ( rho * c_0 * V_wall * K ) / ( 2. * pi )
if optElement[0] == "-x":
B_theta_max = float( optElement[1] ) * pi/180.
V_r_y = V_r * sin( B_theta_max )
V_b_avg = V_r_y + V_b
C_fan = ( rho * c_0 * V_b_avg * K ) / ( 2. * pi )
if optElement[0] == "-z":
AxisOnly = True
PlotNum = 0
if DEBUG:
print "-p", PlotNum
print "-o", OutputFile
print "-f", Frequency
print "-F", Frequency_fan
print "-a", R_i
print "-b", R_o
print "-l", dL
print "-R", F_r
print "-r", range_inc
print "-m", Monopole
print "-v", FlexBaffle
print "-z", AxisOnly
print "-W", delta_wall
print "-x", B_theta_max
return
##----------------------------------------------------------------------------
def OpenFile(file, mode):
try:
# open() will raise an IOError exception if it fails
print "Opening file:", file
fd = open(file,mode)
except IOError, error:
print "OpenFile(): I/O error:", error
sys.exit("OpenFile()")
except OSError, error:
print "OpenFile(): OS error:", error
sys.exit("OpenFile()")
except:
print "OpenFile(): Unexpected error:", sys.exc_info()[0]
raise
sys.exit("OpenFile()")
else: # executed after the try if no exceptions
return fd
##----------------------------------------------------------------------------
## Provide for cmd line invocation
if __name__ == "__main__":
main()