##----------------------------------------------------------------------------
## Name: surf_is.py - create animations and plots
## Purpose: Compute the JONSWAP Wave Spectra, then compute a breaking
## wave height spectrum, a bubble cloud size spectrum, and
## finally an acoustic radiation spectrum from the initial
## bubble oscillation and collective bubble oscillations
## Author: J Park
##
## Created: 11/23/06
##
## Modified:
##----------------------------------------------------------------------------
import os, sys, string, getopt, time
# 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.
# the only reason to use numpy/scipy is for matplotlib graphics
try:
from numpy import * # www.numpy.org numpy.scipy.org
except ImportError:
print "Failed to import numpy."
try:
from pylab import * # matplotlib.sourceforge.net
except ImportError:
print "Failed to import pylab for matplotlib."
from scipy.interpolate.fitpack2 import UnivariateSpline
Version = sys.argv[0]+" : "+str( time.asctime( ( time.gmtime(time.time()) ) ) )
DEBUG = False
##----------------------------------------------------------------------------
# Global variables
g = 9.82
SWH = 0.
peakWaveFreq = 0.
# Default Parameters
plotPolihaleRawData = False
writeOutput = False # override with -z
trainANNData = False # override with -k
runANNData = False # override with -r
createAnimation = False # override with -m
useSphericalCloud = False # override with -a
useInitialBubble = True # override with -j
additionalSpectrum = False # override with -g
windSpeedIncrement = 1. # override with -e
N = 1200 # points in JONSWAP spectrum -n
omegaIncrement = 2.25/N # radian frequency increment of N point wave spectrum
# 2.25 gives a band of 0.25 - 2.5 rad/s (0.04 - 0.4 Hz)
isFreqIncrement = 0.5 # bin width in Hz of IS summed spectrum -q
SpectrumFile = "spect.dat"
InfrasoundFile = "is_surf.dat"
windspeed_knots = 18. # -w input from command line in knots
windspeed = 18. # converted from knots to m/s
fetch_in = 300. # -f input from command line in naut miles
fetch = 300. # converted from naut mile to meters
voidFraction = 0.35 # -v
initialBubbleFreqFactor = 1.1 # frequency multiplier for initial pure air bubble -i
bubbleCylinderFactor = 0.38 # semicylindrical plume radius = fraction wave height -b
bubbleSphereFactor = 0.38 # spherical plume radius = fraction wave height -c
# Global Array's, created with numpy array
JONSWAP_RMK = array([])
waveFreq = array([])
wavelength = array([])
breakHeight = array([])
cylindricalPlumeRadius = array([])
sphericalPlumeRadius = array([])
initialBubbleFreq = array([])
cylinderCloudFreq = array([])
sphereCloudFreq = array([])
isFreq = array([])
isInitBubbleAmp = array([])
isCylinderAmp = array([])
isSphereAmp = array([])
infrasound = array([])
infrasound_dB = array([])
cloudSound_p1T = array([])
cloudSound_p2T = array([])
cloudSound_p3T = array([])
cloudSound_p4T = array([])
cloudSound_p5T = array([])
# global list for ANNTrainData
ANNTrainData = []
ANNRunData = []
# global plotting variables
plotFiles = []
plotSeries = 0
nPlot = 0
# Polihale data
ph_freq = array( [ 2, 4 , 6, 8, 10, 12, 14, 16, 18 ] )
ph_1_68 = array( [ 5, 7.3, 6.8, 5.4, 4, 3, 2.2, 1.5, 1.1 ] )
ph_1_82 = array( [ 7, 9, 7.8, 5.7, 4.2, 3, 2.5, 1.7, 0.9 ] )
ph_1_66 = array( [ 0.1, 1.2, 2.3, 1.9, 1.4, 1.2, 1, 0.7, 0.4 ] )
ph_2_14 = array( [ 7.6, 9.8, 8.8, 6.4, 5.3, 4.3, 3.8, 3, 2.2 ] )
ph_2_09 = array( [ 5.3, 8.3, 7.8, 5.7, 4.5, 3.5, 2.6, 1.8, 1.3 ] )
ph_2_00 = array( [ 1.5, 3.4, 4, 3, 2.3, 1.7, 1.7, 1.4, 1.1 ] )
ph_2_59 = array( [ 7.5, 9.3, 9.1, 6.8, 5.5, 4.8, 4.5, 3.7, 2.9 ] )
ph_2_97 = array( [ 8.8, 10.5, 9.5, 7.3, 5.5, 3.9, 3., 2.2, 1.3 ] )
ph_3_74 = array( [ 15.5, 15.9, 13.9, 11.4, 10.3, 9.2, 8.6, 8, 7.3 ] )
ph_3_83 = array( [ 16, 15.7, 14, 11, 9.5, 8.4, 8, 6.9, 6.2 ] )
ph_4_61 = array( [ 18.1, 17.8, 15.3, 12.8, 11.5, 10.7, 10, 9.4, 8.4 ] )
ph_4_65 = array( [ 17.3, 16.8, 14.2, 11.8, 10.3, 9.4, 8.8, 7.9, 7 ] )
ph_5_33 = array( [ 23, 20.8, 18.1, 15.6, 14.2, 13.2, 12.5, 11.5, 10.5] )
ph_6_15 = array( [ 21, 19, 16.4, 13.8, 12.2, 11.2, 10.5, 9.4, 8.4 ] )
# ordinate arrays for smoothed (spline) Polihale data
y_1_68 = array([])
y_1_82 = array([])
y_1_66 = array([])
y_2_14 = array([])
y_2_09 = array([])
y_2_00 = array([])
y_2_59 = array([])
y_2_97 = array([])
y_3_74 = array([])
y_3_83 = array([])
y_4_61 = array([])
y_4_65 = array([])
y_5_33 = array([])
y_6_15 = array([])
# bubble plume radius increase data: Lamarre & Melville, JASA 95(3)
# the radius at t/T would be R ~ sqrt(A*)/2, with relative values:
# t/T = [0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1]
# A* = [4, 6, 9, 11, 12, 12, 12, 11.5, 11, 10.5, 10]
# R = [1, 1.22, 1.5, 1.66, 1.73, 1.73, 1.73, 1.7, 1.66, 1.62, 1.58]
plumeRadiusIncrease_t_T = array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 ])
plumeRadiusIncrease = array([1.0, 1.22, 1.5, 1.66, 1.73, 1.73, 1.73, 1.7, 1.66, 1.62, 1.58])
plumeRadiusSpline = UnivariateSpline( plumeRadiusIncrease_t_T, plumeRadiusIncrease )
##----------------------------------------------------------------------------
## Main module
def main():
print Version
ParseCmdLine()
global windspeed_knots
SplineData()
if createAnimation:
os.system("rm -f *.png")
# spectrums from 1 to 40 knots, increasing
for wspeed in arange (1., 41., windSpeedIncrement):
windspeed_knots = wspeed
GenerateSpectrums( "a" )
# clear all the data arrays for the next plot
ResetGlobalArrrays()
# spectrums from 40 to 1 knot, decreasing
CopySpectrums( "b" )
# use mplayer mencoder to create an animation in mpg format
os.system("E:\mplayer\mencoder "
"mf://*.png -mf type=png:fps=3 -ovc "
"lavc -lavcopts vcodec=wmv2 -oac copy -o output.mpg")
# delete the plot files
for fname in plotFiles: os.remove(fname)
elif trainANNData or runANNData:
# spectrums from 1 to 40 knots, increasing
for wspeed in arange (1., 41., windSpeedIncrement):
windspeed_knots = wspeed
GenerateSpectrums( None )
UpdateANNData()
# clear all the data arrays for the next iteration
ResetGlobalArrrays()
WriteANNFile()
else:
# single spectrum
GenerateSpectrums( "c" )
show()
if writeOutput:
WriteOutputFiles()
print "Normal Exit"
## Main module
##----------------------------------------------------------------------------
##----------------------------------------------------------------------------
def GenerateSpectrums( figLabelHeader ):
global cloudSound_p1T, cloudSound_p2T, cloudSound_p3T, cloudSound_p4T, cloudSound_p5T
global voidFraction
WaveSpectrum() # compute: JONSWAP_RMK waveFreq wavelength
WaveBreakSpectrum() # compute: breakHeight
WaveBubbleSpectrum() # compute: cylindricalPlumeRadius sphericalPlumeRadius
AcousticBubbleSpectrum() # compute: cylinderCloudFreq sphereCloudFreq
if additionalSpectrum:
cloudSound_p1T = AdditionalBubbleSpectrum( 0.1 )
cloudSound_p2T = AdditionalBubbleSpectrum( 0.2 )
cloudSound_p3T = AdditionalBubbleSpectrum( 0.3 )
cloudSound_p4T = AdditionalBubbleSpectrum( 0.4 )
cloudSound_p5T = AdditionalBubbleSpectrum( 0.5 )
SumAcousticSpectrum() # compute: isFreq infrasound
if figLabelHeader:
CreateGraphics( figLabelHeader )
##----------------------------------------------------------------------------
def CreateGraphics( figLabelHeader ):
global plotFiles, plotSeries, nPlot
# create graphics using pylab (matplotlib)
#-----------------------------------------------------------------
# ocean wave spectra & acoustic radiation freq vs. ocean wave freq
clf()
subplot(2,1,1)
#figure(1)
plot( waveFreq, sqrt( JONSWAP_RMK ),
waveFreq, breakHeight,
waveFreq, initialBubbleFreq,
waveFreq, cylinderCloudFreq, linewidth=2 )
title ( "Wave Spectrum & Acoustic Radiation" )
ylabel( r'$\rm{ S(m) : H_b(m) : f_a(Hz) }$' ) # r - 'raw' string for LaTex
xlabel( "freq (Hz)" )
axis( [ 0, 0.3, 0, 20 ] )
windSpeedLabel = "Wind %.1f knots" % windspeed_knots
SWHLabel = "SWH %.1f m" % SWH
fetchLabel = "Fetch " + str( int(fetch_in) ) + " nm"
figtext( 0.65, 0.75, windSpeedLabel )
figtext( 0.65, 0.72, SWHLabel )
figtext( 0.65, 0.69, fetchLabel )
legend( ("Wave Spectrum","Break Height","Initial Bubble","Collective Bubble") )
#-----------------------------------------------------------------
# acoustic amplitude vs. frequency
subplot(2,1,2)
# figure(2)
# clf()
# line plot of linear summed spectrum
# plot( isFreq, infrasound,
# isFreq, isInitBubbleAmp,
# isFreq, isCylinderAmp, linewidth=2 )
# bar chart (histogram) of linear summed spectrum
# b1 = bar( isFreq, infrasound, width=0.2 )
# bar chart (histogram) of summed spectrum in dB
b1 = bar( isFreq, infrasound_dB, width=isFreqIncrement )
# line plots of polihale data
# If you make multiple lines with one plot command, the kwargs apply
# to all those lines, i.e. linewidth=2, color="r"
yAxisMax = max(infrasound_dB) + 3
if SWH >= 6. and SWH < 10. :
yAxisMax = max(yAxisMax, y_6_15.max())
plot( x_spline, y_6_15, linewidth=2, color="r" )
if plotPolihaleRawData:
plot( ph_freq, ph_6_15, linewidth=2 )
elif SWH >= 5. and SWH < 6. :
yAxisMax = max(yAxisMax, y_5_33.max())
plot( x_spline, y_5_33, linewidth=2, color="r" )
if plotPolihaleRawData:
plot( ph_freq, ph_5_33, linewidth=2 )
elif SWH >= 4. and SWH < 5. :
yAxisMax = max(yAxisMax, y_4_61.max(), y_4_65.max())
plot( x_spline, y_4_61,
x_spline, y_4_65, linewidth=2, color="r" )
if plotPolihaleRawData:
plot( ph_freq, ph_4_61,
ph_freq, ph_4_65, linewidth=2 )
elif SWH >= 3.3 and SWH < 4. :
yAxisMax = max(yAxisMax, y_3_74.max(), y_3_83.max())
plot( x_spline, y_3_74,
x_spline, y_3_83, linewidth=2, color="r" )
if plotPolihaleRawData:
plot( ph_freq, ph_3_74,
ph_freq, ph_3_83, linewidth=2 )
elif SWH >= 2.5 and SWH < 3.3 :
yAxisMax = max(yAxisMax, y_2_59.max(), y_2_97.max())
plot( x_spline, y_2_59,
x_spline, y_2_97, linewidth=2, color="r" )
if plotPolihaleRawData:
plot( ph_freq, ph_2_59,
ph_freq, ph_2_97, linewidth=2 )
elif SWH >= 2.0 and SWH < 2.5 :
yAxisMax = max(yAxisMax, y_2_00.max(), y_2_09.max())
plot( x_spline, y_2_00,
x_spline, y_2_09,
x_spline, y_2_14, linewidth=2, color="r" )
if plotPolihaleRawData:
plot( ph_freq, ph_2_00,
ph_freq, ph_2_09,
ph_freq, ph_2_14, linewidth=2 )
elif SWH >= 1.5 and SWH < 2.0 :
yAxisMax = max(yAxisMax, y_1_66.max(), y_1_68.max())
plot( x_spline, y_1_66,
x_spline, y_1_68,
x_spline, y_1_82, linewidth=2, color="r" )
if plotPolihaleRawData:
plot( ph_freq, ph_1_66,
ph_freq, ph_1_68,
ph_freq, ph_1_82, linewidth=2 )
xlabel( "freq (Hz)" )
ylabel( "amplitude (dB)" )
if SWH >= 2.5 :
legend( ("Polihale","Polihale","Model") )
elif SWH >= 1.5:
legend( ("Polihale","Polihale","Polihale","Model") )
if additionalSpectrum :
figtext( 0.45, 0.42, "Plume Evolution" )
if createAnimation:
yAxisMax = 20
if SWH > 3.:
yAxisMax = 30
axis( [ 0, 20, 0, yAxisMax ] )
else:
axis( [ 0, 20, max( 0, infrasound_dB[-1] - 5 ), yAxisMax ] )
xticks( [0,2,4,6,8,10,12,14,16,18,20] )
#-----------------------------------------------------------------
# create file name that will be read in order created when
# mencoder is called to create the animation
nPlot += 1
if nPlot > 9 and nPlot < 20:
plotSeries = 1
elif nPlot > 19 and nPlot < 30:
plotSeries = 2
elif nPlot > 29 and nPlot < 40:
plotSeries = 3
elif nPlot > 39 and nPlot < 50:
plotSeries = 4
elif nPlot > 49 and nPlot < 60:
plotSeries = 5
elif nPlot > 59 and nPlot < 70:
plotSeries = 6
elif nPlot > 69 and nPlot < 80:
plotSeries = 7
elif nPlot > 79 and nPlot < 90:
plotSeries = 8
elif nPlot > 89 and nPlot < 100:
plotSeries = 9
plotName = "fig" + figLabelHeader + str(plotSeries) + "_" + str(nPlot) + ".png"
plotFiles.append( plotName )
# savefig will output to these formats:
# eps, jpeg, pdf, png, ps, svg
# the format is deduced by the filename extension (.png is default)
savefig( plotName )
##----------------------------------------------------------------------------
def CopySpectrums( figLabelHeader ):
global plotFiles, plotSeries, nPlot
plotSeries = 0
nPlot = 0
plotFiles2 = []
for i in range ( len(plotFiles)-1, -1, -1 ) :
nPlot += 1
if nPlot > 9 and nPlot < 20:
plotSeries = 1
elif nPlot > 19 and nPlot < 30:
plotSeries = 2
elif nPlot > 29 and nPlot < 40:
plotSeries = 3
elif nPlot > 39 and nPlot < 50:
plotSeries = 4
elif nPlot > 49 and nPlot < 60:
plotSeries = 5
elif nPlot > 59 and nPlot < 70:
plotSeries = 6
elif nPlot > 69 and nPlot < 80:
plotSeries = 7
elif nPlot > 79 and nPlot < 90:
plotSeries = 8
elif nPlot > 89 and nPlot < 100:
plotSeries = 9
plotName = "fig" + figLabelHeader + str(plotSeries) + "_" + str(nPlot) + ".png"
plotFiles2.append( plotName )
os.system("cp %s %s" % (plotFiles[i], plotName) )
for i in range ( len(plotFiles2) ) :
plotFiles.append( plotFiles2[i] )
if DEBUG:
print plotFiles
##----------------------------------------------------------------------------
def SplineData():
global x_spline, y_1_68, y_1_82, y_1_66, y_2_14, y_2_09, y_2_00
global y_2_59, y_2_97, y_3_74, y_3_83, y_4_61, y_4_65, y_5_33, y_6_15
# spline smoothing for polihale data
x_spline = linspace(2, 18, 200)
spline = UnivariateSpline( ph_freq, ph_1_68, s=0.01 )
y_1_68 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_1_82, s=0.01 )
y_1_82 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_1_66, s=0.01 )
y_1_66 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_2_14, s=0.01 )
y_2_14 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_2_09, s=0.01 )
y_2_09 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_2_00, s=0.01 )
y_2_00 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_2_59, s=0.01 )
y_2_59 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_2_97, s=0.01 )
y_2_97 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_3_74, s=0.01 )
y_3_74 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_3_83, s=0.01 )
y_3_83 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_4_61, s=0.01 )
y_4_61 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_4_65, s=0.01 )
y_4_65 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_5_33, s=0.01 )
y_5_33 = spline( x_spline )
spline = UnivariateSpline( ph_freq, ph_6_15, s=0.01 )
y_6_15 = spline( x_spline )
##----------------------------------------------------------------------------
def SumAcousticSpectrum():
global isFreq, isSphereAmp, isCylinderAmp, isInitBubbleAmp
global infrasound, infrasound_dB
# this is the big unknown
# how to assign spectral amplitudes to the acoustic radiation frequencies?
# setup an acoustic radiation frequency 'x-axis'
freqMin = 0.001
freqMax = 20
numFreq = int( ceil ( (freqMax - freqMin) / isFreqIncrement ) )
freq = freqMin
# create zero'd arrays of the right length for accumulation
for i in range( numFreq ):
isFreq = append( isFreq, freq )
infrasound = append( infrasound, 0 )
infrasound_dB = append( infrasound_dB, 0 )
isSphereAmp = append( isSphereAmp, 0 )
isCylinderAmp = append( isCylinderAmp, 0 )
isInitBubbleAmp = append( isInitBubbleAmp, 0 )
freq = freq + isFreqIncrement
# 'integrate' amplitude contributions of acoustic radiation frequency
# across the ocean wave frequency domain of single bubble and cloud
# frequencies, sum both to represent the spectral shape of the infrasound
# assume a 1/f amplitide scaling of the radiation frequencies, and
# modulate this 1/f by a coefficient A
# if we were adding 1 instead of A * ( 1 / isFreq[j] ) that would
# simply bin (create histogram) of single bubble and cloud frequencies
A = 2.
for i in range(N):
for j in range( numFreq ):
if sphereCloudFreq[i] > (isFreq[j] - isFreqIncrement/2) and \
sphereCloudFreq[i] < (isFreq[j] + isFreqIncrement/2) :
isSphereAmp[j] = isSphereAmp[j] + A / isFreq[j]
if useSphericalCloud:
infrasound[j] = infrasound[j] + A / isFreq[j]
break
for j in range( numFreq ):
if cylinderCloudFreq[i] > (isFreq[j] - isFreqIncrement/2) and \
cylinderCloudFreq[i] < (isFreq[j] + isFreqIncrement/2) :
isCylinderAmp[j] = isCylinderAmp[j] + A / isFreq[j]
if not useSphericalCloud:
infrasound[j] = infrasound[j] + A / isFreq[j]
break
if useInitialBubble:
for j in range( numFreq ):
if initialBubbleFreq[i] > (isFreq[j] - isFreqIncrement/2) and \
initialBubbleFreq[i] < (isFreq[j] + isFreqIncrement/2) :
isInitBubbleAmp[j] = isInitBubbleAmp[j] + A / isFreq[j]
infrasound[j] = infrasound[j] + A / isFreq[j]
break
if additionalSpectrum:
for j in range( numFreq ):
if cloudSound_p1T[i] > (isFreq[j] - isFreqIncrement/2) and \
cloudSound_p1T[i] < (isFreq[j] + isFreqIncrement/2) :
infrasound[j] = infrasound[j] + A / isFreq[j]
break
for j in range( numFreq ):
if cloudSound_p2T[i] > (isFreq[j] - isFreqIncrement/2) and \
cloudSound_p2T[i] < (isFreq[j] + isFreqIncrement/2) :
infrasound[j] = infrasound[j] + A / isFreq[j]
break
for j in range( numFreq ):
if cloudSound_p3T[i] > (isFreq[j] - isFreqIncrement/2) and \
cloudSound_p3T[i] < (isFreq[j] + isFreqIncrement/2) :
infrasound[j] = infrasound[j] + A / isFreq[j]
break
for j in range( numFreq ):
if cloudSound_p4T[i] > (isFreq[j] - isFreqIncrement/2) and \
cloudSound_p4T[i] < (isFreq[j] + isFreqIncrement/2) :
infrasound[j] = infrasound[j] + A / isFreq[j]
break
for j in range( numFreq ):
if cloudSound_p5T[i] > (isFreq[j] - isFreqIncrement/2) and \
cloudSound_p5T[i] < (isFreq[j] + isFreqIncrement/2) :
infrasound[j] = infrasound[j] + A / isFreq[j]
break
# convert summed spectrum to log scale (dB)
for i in range( numFreq ):
if infrasound[i] < 1E-4:
infrasound_dB[i] = 0
else:
infrasound_dB[i] = 10. * log10( infrasound[i] )
# since the amplitudes are dependent on N, normalize to polihale data
# use the 10 Hz value of the polihale data
i10Hz = ph_freq.searchsorted( 10 )
if i10Hz >= len( ph_freq ) :
print "Error: SumAcousticSpectrum(): invalid index"
sys.exit(1)
polihaleAmp10Hz = -1. # SWH < 0.8
if additionalSpectrum :
polihaleAmp10Hz = -7. # SWH < 0.8
if SWH >= 6. and SWH < 10. :
polihaleAmp10Hz = ph_6_15[ i10Hz ]
elif SWH >= 5. and SWH < 6. :
polihaleAmp10Hz = ph_5_33[ i10Hz ]
elif SWH >= 4. and SWH < 5. :
polihaleAmp10Hz = max( ph_4_61[ i10Hz ], ph_4_65[ i10Hz ] )
elif SWH >= 3.3 and SWH < 4. :
polihaleAmp10Hz = max( ph_3_74[ i10Hz ], ph_3_83[ i10Hz ] )
elif SWH >= 2.5 and SWH < 3.3 :
polihaleAmp10Hz = max( ph_2_59[ i10Hz ], ph_2_97[ i10Hz ] )
elif SWH >= 2.0 and SWH < 2.5 :
polihaleAmp10Hz = max( ph_2_00[ i10Hz ], ph_2_09[ i10Hz ] )
elif SWH >= 1.5 and SWH < 2.0 :
polihaleAmp10Hz = max( ph_1_66[ i10Hz ], ph_1_68[ i10Hz ] )
elif SWH >= 0.8 and SWH < 1.5 :
polihaleAmp10Hz = max( ph_1_66[ i10Hz ], ph_1_68[ i10Hz ] )
# difference between summed spectrum and polihale @ 10 Hz
i10Hz = isFreq.searchsorted( 10 )
delta_dB = polihaleAmp10Hz - infrasound_dB[i10Hz]
#print "delta_dB = ", polihaleAmp10Hz, " - ", infrasound_dB[i10Hz], " = ", delta_dB
for i in range( numFreq ):
if infrasound_dB[i] > 0.1:
infrasound_dB[i] = infrasound_dB[i] + delta_dB
##----------------------------------------------------------------------------
def AdditionalBubbleSpectrum( t_T ):
# create additional contributions based on the reduction of void
# fraction and increase of plume cross section as time progresses
# input is t/T = fraction of wave period
# create a new array returned by this function
arrayIn = array([])
for i in range( len( cylindricalPlumeRadius ) ):
if cylindricalPlumeRadius[i] < 0.1 :
arrayIn = append( arrayIn, 999 )
continue
# based on Lamarre & Melville, JASA 95(3) p1317, figure 6(c)
# the void gradient is dv/dt = -0.56 for t/T < 0.5
reduced_voidFraction = voidFraction - 0.56 * t_T
if reduced_voidFraction < 0.01:
reduced_voidFraction = 0.01
# based on Lamarre & Melville, JASA 95(3) p1317, figure 6(b)
dx = plumeRadiusIncrease_t_T[1] - plumeRadiusIncrease_t_T[0]
interpolatePoint = array([t_T])
r = plumeRadiusSpline( interpolatePoint )
increasedPlumeRadius = cylindricalPlumeRadius[i] * r
# Breathing-mode oscillation of semicylindrical bubble plume of
# radius bubble[i], from: Lowen & Melville, JASA 95(3), p1329
f_R = sqrt( 101325 / ( 1000 * reduced_voidFraction ) ) \
/ ( 2.5 * increasedPlumeRadius )
arrayIn = append( arrayIn, f_R )
return arrayIn
##----------------------------------------------------------------------------
def AcousticBubbleSpectrum():
global cylinderCloudFreq, sphereCloudFreq, initialBubbleFreq
for i in range( len( sphericalPlumeRadius ) ):
if sphericalPlumeRadius[i] < 0.1 :
sphereCloudFreq = append( sphereCloudFreq, 999 )
cylinderCloudFreq = append( cylinderCloudFreq, 999 )
initialBubbleFreq = append( initialBubbleFreq, 999 )
continue
# Lowest Eigenfrequency of collective bubble oscillation for cloud
# of radius sphericalPlumeRadius[i], from: Carey 1993, reported in Deane
# JASA 102(5), 1997
x = pow( 3 * 101325 / ( 1000 * voidFraction ), 0.5 )
f_0 = (1 / (2 * pi * sphericalPlumeRadius[i])) * x
sphereCloudFreq = append( sphereCloudFreq, f_0 )
# Breathing-mode oscillation of semicylindrical bubble plume of
# radius bubble[i], from: Lowen & Melville, JASA 95(3), p1329
f_R = sqrt( 101325 / ( 1000 * voidFraction ) ) / ( 2.5 * cylindricalPlumeRadius[i] )
cylinderCloudFreq = append( cylinderCloudFreq, f_R )
# JP: analytical comparison shows f_R = 1.45 * f_0
# However, the size range for this expression is 0.1 < R_c < 1 (m)
# with void fractions less than 0.1
# So we are at scales clearly outside their recommendation...
# Nonetheless, observations show that the real-world plumes are
# semicylindrical rather than spherical...
# and therefore have a higher frequency than a spherical plume
# will use cylinderCloudFreq as default, user can override with -a
# volumetric pulsation of the cylinder of pure air with a slightly higher
# frequency than the plume... apply a simple scale factor
initialBubbleFreq = append( initialBubbleFreq,
initialBubbleFreqFactor * cylinderCloudFreq[-1] )
##----------------------------------------------------------------------------
def WaveBubbleSpectrum():
global cylindricalPlumeRadius, sphericalPlumeRadius
# convert the breaking wave height spectrum into a spectrum
# of semicylindrical bubble and spherical plume radius based on
# scale factors applied to the breaking wave height
for i in range( len(breakHeight) ):
if breakHeight[i] < 0.1 :
cynBubbleRadius = 0
sphereBubbleRadius = 0
else:
cynBubbleRadius = breakHeight[i] * bubbleCylinderFactor
sphereBubbleRadius = breakHeight[i] * bubbleSphereFactor
cylindricalPlumeRadius = append ( cylindricalPlumeRadius, cynBubbleRadius )
sphericalPlumeRadius = append ( sphericalPlumeRadius, sphereBubbleRadius )
#print bubbleRadius,
##----------------------------------------------------------------------------
def WaveBreakSpectrum():
global breakHeight
# convert the deepwater wave spectrum into a breaking wave height
# spectrum by application of the 'breaking index'
for i in range( len(JONSWAP_RMK) ):
deepwaterHeight = sqrt( JONSWAP_RMK[i] )
WaveSteepness = deepwaterHeight / wavelength[i]
breakIndex = WaveBreakIndex( WaveSteepness )
breakHeight = append( breakHeight, deepwaterHeight * breakIndex )
##----------------------------------------------------------------------------
def WaveBreakIndex( WaveSteepness ):
# Compute an 'average' breaker height index from figure 2-65 of the SPM
if WaveSteepness < 0.002:
return 1.8
if WaveSteepness > 0.1:
return 1.
# upper bound is for beach slope of m = 1/10 = 0.1
hiIndex = 0.578422 * pow( WaveSteepness, -0.239253 )
# lower bound is for beach slope of m = 1/50 = 0.02
loIndex = 0.5403 * pow( WaveSteepness, -0.216176 )
# since we don't have beach slope, take an average
breakIndex = ( hiIndex + loIndex ) / 2
# limit to the plunging breaker range of figure 2-65
if breakIndex < 1.2 : breakIndex = 1.2
elif breakIndex > 1.8 : breakIndex = 1.8
return breakIndex
##----------------------------------------------------------------------------
def WaveSpectrum():
global fetch, SWH, peakWaveFreq
global JONSWAP_RMK, waveFreq, wavelength
omega = 0.1 # initial frequency (rad/s)
windspeed = windspeed_knots / 1.9425 # convert knots to m/s
fetch = fetch_in * 1853 # convert nautical miles to m
print "fetch=%6.0f (m) %4.0f (nm) windspeed=%5.2f (m/s) %5.2f (knots)" % \
(fetch, fetch/1853, windspeed, windspeed_knots)
# Parameters for JONSWAP_RMK spectrum from 'A Probalistic Description of Deepwater
# Whitecaps', Robert M. Kennedy, PhD Dissertaion, Nova University, May 1978, p 5-8
nd_fetch = g * fetch / pow( windspeed, 2 ) # non-dimensional fetch
Omega = 22. * (g / windspeed) * pow( nd_fetch, -0.33 ) # max amplitude freq
peakWaveFreq = Omega / (2*pi)
epsilon = 0.076 * pow( nd_fetch , -0.22 ) # Phillips constant (alpha)
lambda_peak = 2*pi * g / pow( Omega, 2 ) # max amplitude wavelength
#print "peakWaveFreq=%6.5f nd_fetch=%6.0f epsilon=%7.6f lambda_peak=%4.0f" % \
# ( peakWaveFreq, nd_fetch, epsilon, lambda_peak )
for i in range(N):
f = omega/(2*pi)
# Compute the JONSWAP_RMK spectrum for each frequency omega
S1 = epsilon * pow( g, 2 ) * pow( omega, -5 )
S2 = exp( -(5/4) * pow( omega/Omega, -4) )
S3 = pow( 3.3, exp( -(pow(omega - Omega, 2)/(0.0128 * pow(Omega, 2))) ) )
S_RMK = S1 * S2 * S3
# save the spectral amplitudes in the arrays
JONSWAP_RMK = append( JONSWAP_RMK, S_RMK )
# save the wave frequency in Hz
waveFreq = append( waveFreq, f )
# save the wavelength
wavelength = append( wavelength, ( g * pow(1/f,2) )/( 2*pi ) )
# increment the frequency
omega = omega + omegaIncrement
if DEBUG:
print "waveFreq: ", waveFreq
print "wavelength: ", wavelength
print "JONSWAP_RMK: ", JONSWAP_RMK
# Compute significant wave height
SWH = 0.
for i in range(N):
SWH = SWH + JONSWAP_RMK[i] * omegaIncrement
SWH = SWH / (2*pi)
SWH = 4 * sqrt( SWH )
print "Significant wave height=%5.2f (m) %5.2f (ft)" % (SWH, SWH/.3041)
##----------------------------------------------------------------------------
def ResetGlobalArrrays():
global JONSWAP_RMK, waveFreq, wavelength, breakHeight
global cylindricalPlumeRadius, sphericalPlumeRadius, initialBubbleFreq
global cylinderCloudFreq, sphereCloudFreq, isFreq
global isInitBubbleAmp, isCylinderAmp
global isSphereAmp, infrasound, infrasound_dB
JONSWAP_RMK = array([])
waveFreq = array([])
wavelength = array([])
breakHeight = array([])
cylindricalPlumeRadius = array([])
sphericalPlumeRadius = array([])
initialBubbleFreq = array([])
cylinderCloudFreq = array([])
sphereCloudFreq = array([])
isFreq = array([])
isInitBubbleAmp = array([])
isCylinderAmp = array([])
isSphereAmp = array([])
infrasound = array([])
infrasound_dB = array([])
##----------------------------------------------------------------------------
def WriteOutputFiles():
##
## Open the output files
if SpectrumFile and InfrasoundFile:
fdSpectrum = OpenFile(SpectrumFile, 'w')
fdInfrasound = OpenFile(InfrasoundFile, 'w')
else:
print "Missing output file specification.\n"
usage()
return
# Write the spectrums to the output file
# Note the wave spectrum is reduced to units of (m) via sqrt()
fdSpectrum.write("WaveFrequency(Hz) JONSWAP(m) breakHeight(m) "
"CylinderFrequency(Hz) InitBubbleFrequency(Hz) "
"PeakFrequency=%7.6f(Hz) SWH=%3.1f(m) void=%3.2f\n" % \
(peakWaveFreq, SWH, voidFraction))
for i in range(N):
# cylindricalK = pow( (2 * pi * cylinderCloudFreq[i]), 2 ) / g
# oceanWaveNumber = pow( (2 * pi * waveFreq[i]), 2 ) / g
fdSpectrum.write("%.5f %.3f %.3f %.3f %.3f\n" % \
(waveFreq[i], sqrt(JONSWAP_RMK[i]), breakHeight[i], \
cylinderCloudFreq[i], initialBubbleFreq[i] ))
# Write the infrasound to the output file
fdInfrasound.write("Frequency(Hz) IS(dB) IS_Sphere IS_Cylinder IS_InitBubble "
"PeakFrequency=%7.6f(Hz) SWH=%3.1f(m) void=%3.2f\n" % \
(peakWaveFreq, SWH, voidFraction))
for i in range( len(infrasound) ):
fdInfrasound.write("%.4f %.2f %.0f %.0f %.0f\n" % \
( isFreq[i], infrasound_dB[i], isSphereAmp[i], \
isCylinderAmp[i], isInitBubbleAmp[i] ))
print "Wrote %d lines to %s" % (len(JONSWAP_RMK), SpectrumFile)
print "Wrote %d lines to %s" % (len(infrasound_dB), InfrasoundFile)
## Close the files
fdSpectrum.close()
fdInfrasound.close()
##----------------------------------------------------------------------------
def UpdateANNData():
global ANNTrainData, ANNRunData
# ANNTrainData is a list of strings representing the input/output
# pairs of infrasound_dB (B(wa)) and the breaking wave spectrum
# breakHeight (Hb(w)). Each neuron will represent a frequency,
# so if there are 20 values of infrasound_dB, then each value
# will correspond to a specific acoustic frequency. The actual
# frequency is not needed by the ANN, just the spectral amplitudes.
# Two strings are added on each call, one Input, one Output.
# The whole list is dumped to a file in WriteANNFiles()
# down sample the infrasound_DB spectrum into
# 20 points evenly spaced from 1-19 Hz with a cubic spline
omega_output = linspace(1., 19., 20)
spline = UnivariateSpline( isFreq, infrasound_dB )
IS_Spectrum = spline( omega_output )
# Create a string for the input data
inputData = "I:w%2d:SWH%3.1f: " % (windspeed_knots, SWH)
for i in range( len(IS_Spectrum) ):
# take abs to avoid values such as -0.01 dB
inputData += "%.2f " % abs( IS_Spectrum[i] )
inputData += "\n"
if trainANNData:
# down sample the breakHeight spectrum from N points into
# 20 points evenly spaced from 0.035 - 0.2 Hz with a cubic spline
omega_output = linspace(0.035, 0.2, 20)
spline = UnivariateSpline( waveFreq, breakHeight, s=0.05 )
breakSpectrum = spline( omega_output )
# Create a string for the output data
outputData = "O:w%2d:SWH%3.1f: " % (windspeed_knots, SWH)
for i in range( len( breakSpectrum ) ):
# take abs to avoid values such as -0.01 dB
outputData += "%.3f " % abs( breakSpectrum[i] )
outputData += "\n"
# Append this line to both the input and output
ANNTrainData.append(inputData);
ANNTrainData.append(outputData);
elif runANNData:
# Append only to the input, there is no output needed
ANNRunData.append(inputData);
else:
print "Error: UpdateANNData() invalid state"
sys.exit(1)
##----------------------------------------------------------------------------
def WriteANNFile():
##
if trainANNData:
ANNFile = "surf_is_TrainANN.dat"
ANNData = ANNTrainData
elif runANNData:
ANNFile = "surf_is_RunANN.dat"
ANNData = ANNRunData
## Open the output file
fdANN = OpenFile(ANNFile, 'w')
for i in range( len(ANNData) ):
fdANN.write( ANNData[i] )
print "Wrote %d lines to %s" % (len(ANNData), ANNFile)
## Close the file
fdANN.close()
##----------------------------------------------------------------------------
def usage():
print "Usage: ", sys.argv[0]
print "\t -h (help)"
print "\t -d (debug)"
print "\t -z (write output)"
print "\t -w windspeed (knots)"
print "\t -f fetch (nautmiles)"
print "\t -o spectrumOut"
print "\t -s infrasoundOut"
print "\t -n Nspectrum"
print "\t -b bubbleCylinderFactor"
print "\t -c bubbleSphereFactor"
print "\t -i initialBubbleFreqFactor"
print "\t -v voidFraction"
print "\t -a (use spherical bubble)"
print "\t -j (don't use initial bubble)"
print "\t -q isFreqIncrement"
print "\t -e windSpeedIncrement"
print "\t -g (additional Spectrum)"
print "\t -k (write ANN train data)"
print "\t -r (write ANN run data)"
print "\t -m (create animation)"
##----------------------------------------------------------------------------
def ParseCmdLine():
##
global N, useSphericalCloud, isFreqIncrement, omegaIncrement
global SpectrumFile, InfrasoundFile, windspeed_knots, fetch_in, voidFraction
global bubbleCylinderFactor, bubbleSphereFactor, initialBubbleFreqFactor
global DEBUG, writeOutput, windSpeedIncrement, additionalSpectrum
global useInitialBubble, trainANNData, runANNData, createAnimation
try:
## Parse the command line arguments in sys.argv via getopt()
##
## getopt (args, options[, long_options])
## args is the argument list to be parsed,
## options is the string of option letters that the script wants to
## recognize, with options that require an argument followed by a colon.
## long_options, if specified, must be a list of strings with the names
## of the long options which should be supported.
## The return value consists of two lists: the first is a list of
## (option, value) tuple pairs; the second is the list of program
## arguments left after the option list was stripped.
##
opts, args = getopt.getopt(sys.argv[1:], "hdzw:f:o:s:n:b:c:i:v:ajq:e:gkmr")
except getopt.GetoptError, error:
print "ParseCmdLine(): GetoptError: ",error
return
if len(opts) == 0:
usage()
return
if DEBUG:
print "Options are:",opts,"Length:",len(opts)
print "Arguments are:",args
if len(opts) > 1:
print "opts[0][0]",opts[0][0]," [0][1]", opts[0][1],
print " [1][0]", opts[1][0], " [1][1]",opts[1][1]
argCount = 0
for optElement in opts:
argCount = argCount + 1
if optElement[0] == "-h":
usage()
sys.exit(1)
if optElement[0] == "-d":
DEBUG = True
if optElement[0] == "-z":
writeOutput = True
if optElement[0] == "-e":
windSpeedIncrement = float(optElement[1])
if optElement[0] == "-g":
additionalSpectrum = True
if optElement[0] == "-k":
trainANNData = True
if optElement[0] == "-r":
runANNData = True
if optElement[0] == "-m":
createAnimation = True
if optElement[0] == "-o":
SpectrumFile = optElement[1]
if optElement[0] == "-s":
InfrasoundFile = optElement[1]
if optElement[0] == "-w":
windspeed_knots = float(optElement[1])
if optElement[0] == "-f":
fetch_in = float(optElement[1])
if optElement[0] == "-b":
bubbleCylinderFactor = float(optElement[1])
if optElement[0] == "-c":
bubbleSphereFactor = float(optElement[1])
if optElement[0] == "-i":
initialBubbleFreqFactor = float(optElement[1])
if optElement[0] == "-v":
voidFraction = float(optElement[1])
if optElement[0] == "-a":
useSphericalCloud = True
if optElement[0] == "-j":
useInitialBubble = False
if optElement[0] == "-n":
N = int(optElement[1])
omegaIncrement = 2.25/N
if optElement[0] == "-q":
isFreqIncrement = float(optElement[1])
if trainANNData and createAnimation:
print "Error: can't set both -k (trainANNData) and " \
"-m (createAnimation)"
sys.exit(1)
if DEBUG:
print "opts element %d:" % (argCount), optElement
print "-o", SpectrumFile, "-s", InfrasoundFile, "-w", windspeed_knots, \
"-f", fetch_in, "-b", bubbleCylinderFactor, "-c", bubbleSphereFactor, \
"-i", initialBubbleFreqFactor, "-v", voidFraction, \
"-a", useSphericalCloud, "-n", N, "-q", isFreqIncrement, \
"-z", writeOutput, "-e", windSpeedIncrement, "-g", additionalSpectrum, \
"-j", useInitialBubble, "-k", trainANNData, "-r", runANNData, \
"-m", createAnimation
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()