Commit f0efe60d authored by David Van Komen's avatar David Van Komen

Did some basic cleanup, realized it was built for windows, proj. suspended

parent 3f477c5d
# KRAKEN - A Normal Mode Program in Python - BYU Fork
**Turns out this library is written for Windows execution.**
This repository is a fork of the Python source code for KRAKEN. KRAKEN is a
normal mode program designed to find normal modes in underwater acoustics.
......@@ -27,4 +29,4 @@ Over subsequent years, KRAKEN was greatly extended, with options for modeling oc
This report addresses the need for a more complete user's guide to supplement the on-line help files. The first chapters give a fairly technical description of the mathematical and numerical basis of the model. Additional chapters give a simpler description of its use and installation in a manner that is accessible to less scientifically-oriented readers.
_Found on [HLS Research](http://oalib.hlsresearch.com/Modes/AcousticsToolbox/manual_html/kraken.html)_
\ No newline at end of file
_Found on [HLS Research](http://oalib.hlsresearch.com/Modes/AcousticsToolbox/manual_html/kraken.html)_
......@@ -3,8 +3,13 @@
# Bellhop3D: ray tracing test
# Mexilhoeira Grande, Dom Dez 25 22:50:36 WET 2016
# Written by Tordar
#
#======================================================================
#
# Modified by David Van Komen (BYU)
#
#=======================================================================
# TODO: Fix these messy imports
# That also means fixing all of the calls to functions...
# commandline: python bellhop3d_wedge_rays.py
# python : execfile('bellhop3d_wedge_rays.py')
......@@ -19,16 +24,17 @@ from plotray3d import *
case_title = '''Wedge problem'''
freq = 122.0
ray_step = 1.0
slope = 4.55; k = tan( slope*pi/180 )
z0 = 44.4
x1 = -z0/k
x2 = -x1
z1 = k*x1 + z0
z2 = k*x2 + z0
xmax = max([x1,x2])
zmax = max([z1,z2])
freq = 122.0
ray_step = 1.0
slope = 4.55
k = tan( slope*pi/180 )
z0 = 44.4
x1 = -z0/k
x2 = -x1
z1 = k*x1 + z0
z2 = k*x2 + z0
xmax = max([x1, x2])
zmax = max([z1, z2])
#==================================================================
#
......@@ -36,17 +42,28 @@ zmax = max([z1,z2])
#
#==================================================================
xs = array([0.0])
ys = array([0.0])
zs = array([8.3])
box = array([ray_step, xmax/1000.0, 5.0, zmax]) # [step(m) xbox(km) ybox(km) zbox(m)]
thetas = array([5,6,-27.0,27.0])
phi = array([3,46,90.0,110.0])
bearings = array([2,0.0,1.0])
p = []
comp = []
source_data = {"xs":xs, "ys":ys, "zs":zs, "box":box, "f":freq, "thetas":thetas, "phi":phi, "bearings":bearings, "p":p, "comp":comp}
xs = array([0.0])
ys = array([0.0])
zs = array([8.3])
box = array([ray_step, xmax/1000.0, 5.0, zmax]) # [step(m) xbox(km) ybox(km) zbox(m)]
thetas = array([5,6,-27.0,27.0])
phi = array([3,46,90.0,110.0])
bearings = array([2,0.0,1.0])
p = []
comp = []
source_data = {
"xs": xs,
"ys": ys,
"zs": zs,
"box": box,
"f": freq,
"thetas": thetas,
"phi": phi,
"bearings": bearings,
"p": p,
"comp": comp
}
#==================================================================
#
......@@ -62,7 +79,14 @@ yati = zeros(1)
zati = zeros(1)
p = zeros(1)
surface_data = {"itype":itype,"xati":xati,"yati":yati,"zati":zati,"p":p,"units":aunits}
surface_data = {
"itype": itype,
"xati": xati,
"yati": yati,
"zati": zati,
"p": p,
"units": aunits
}
#==================================================================
#
......@@ -72,11 +96,15 @@ surface_data = {"itype":itype,"xati":xati,"yati":yati,"zati":zati,"p":p,"units":
c0 = 1488.2
z = array([0.0,zmax])
c = array([c0,c0])
z = array([0.0, zmax])
c = array([c0, c0])
r = zeros(1)
ssp_data = {"r":r,"z":z,"c":c}
ssp_data = {
"r": r,
"z": z,
"c": c
}
#==================================================================
#
......@@ -84,14 +112,21 @@ ssp_data = {"r":r,"z":z,"c":c}
#
#==================================================================
bunits = '''W''' # Bottom absorption units: (dB/m)*kHZ
itype = '''R''' # Not used, but required...
xbty = array([x1,x2])/1000.0
ybty = array([-5.0,5.0])
zbty = array([[z1,z1],[z2,z2]])
p = array([1700.0,0.0,1.99,0.5])
bottom_data = {"itype":itype,"xbty":xbty,"ybty":ybty,"zbty":zbty,"p":p,"units":bunits}
bunits = '''W''' # Bottom absorption units: (dB/m)*kHZ
itype = '''R''' # Not used, but required...
xbty = array([x1, x2])/1000.0
ybty = array([-5.0, 5.0])
zbty = array([[z1, z1], [z2, z2]])
p = array([1700.0, 0.0, 1.99, 0.5])
bottom_data = {
"itype": itype,
"xbty": xbty,
"ybty": ybty,
"zbty": zbty,
"p": p,
"units": bunits
}
#==================================================================
#
......@@ -109,12 +144,20 @@ nr = 11
rd = array([1.0,4.0])
r = array([0.0,5.0])
options = {"options1":options1,"options2":options2,"options3":options3,"options4":options4,\
"nrd":nrd,"nr":nr,"rd":rd,"r":r}
options = {
"options1": options1,
"options2": options2,
"options3": options3,
"options4": options4,
"nrd": nrd,
"nr": nr,
"rd": rd,
"r": r
}
print("Writing environmental file...")
wbellhop3denv('wedge',case_title,source_data,surface_data,ssp_data,bottom_data,options)
wbellhop3denv('wedge', case_title, source_data, surface_data, ssp_data, bottom_data, options)
print( "Running Bellhop3D..." )
......
......@@ -4,7 +4,12 @@
# Mexilhoeira Grande, Dom Dez 25 22:55:31 WET 2016
# Written by Tordar
#
#======================================================================
# Modified by David Van Komen (BYU)
#
#=======================================================================
# TODO: Fix these messy imports
# That also means fixing all of the calls to functions...
# commandline: python bellhop3d_wedge_tl.py
# python : execfile('bellhop3d_wedge_tl.py')
......@@ -19,16 +24,17 @@ from readshd import *
case_title = '''Wedge problem'''
freq = 122.0
ray_step = 1.0
slope = 4.55; k = tan( slope*pi/180 )
z0 = 44.4
x1 = -z0/k
x2 = -x1
z1 = k*x1 + z0
z2 = k*x2 + z0
xmax = max([x1,x2])
zmax = max([z1,z2])
freq = 122.0
ray_step = 1.0
slope = 4.55
k = tan( slope*pi/180 )
z0 = 44.4
x1 = -z0/k
x2 = -x1
z1 = k*x1 + z0
z2 = k*x2 + z0
xmax = max([x1,x2])
zmax = max([z1,z2])
#==================================================================
#
......@@ -36,17 +42,28 @@ zmax = max([z1,z2])
#
#==================================================================
xs = array([0.0])
ys = array([0.0])
zs = array([8.3])
box = array([ray_step, xmax/1000.0, 8.0, zmax]) # [step(m) xbox(km) ybox(km) zbox(m)]
thetas = array([401,6,-24.0,24.0])
phi = array([501,46,90.0,110.0])
bearings = array([1,90.0,0.0])
p = []
comp = []
source_data = {"xs":xs, "ys":ys, "zs":zs, "box":box, "f":freq, "thetas":thetas, "phi":phi, "bearings":bearings, "p":p, "comp":comp}
xs = array([0.0])
ys = array([0.0])
zs = array([8.3])
box = array([ray_step, xmax/1000.0, 8.0, zmax]) # [step(m) xbox(km) ybox(km) zbox(m)]
thetas = array([401, 6, -24.0, 24.0])
phi = array([501, 46, 90.0, 110.0])
bearings = array([1, 90.0, 0.0])
p = []
comp = []
source_data = {
"xs": xs,
"ys": ys,
"zs": zs,
"box": box,
"f": freq,
"thetas": thetas,
"phi": phi,
"bearings": bearings,
"p": p,
"comp": comp
}
#==================================================================
#
......@@ -54,15 +71,22 @@ source_data = {"xs":xs, "ys":ys, "zs":zs, "box":box, "f":freq, "thetas":thetas,
#
#==================================================================
aunits = ''' '''
itype = ''' '''
itype = []
xati = zeros(1)
yati = zeros(1)
zati = zeros(1)
p = zeros(1)
surface_data = {"itype":itype,"xati":xati,"yati":yati,"zati":zati,"p":p,"units":aunits}
aunits = ''' '''
itype = ''' '''
itype = []
xati = zeros(1)
yati = zeros(1)
zati = zeros(1)
p = zeros(1)
surface_data = {
"itype": itype,
"xati": xati,
"yati": yati,
"zati": zati,
"p": p,
"units": aunits
}
#==================================================================
#
......@@ -72,11 +96,15 @@ surface_data = {"itype":itype,"xati":xati,"yati":yati,"zati":zati,"p":p,"units":
c0 = 1488.2
z = array([0.0,zmax])
c = array([c0,c0])
z = array([0.0, zmax])
c = array([c0, c0])
r = zeros(1)
ssp_data = {"r":r,"z":z,"c":c}
ssp_data = {
"r": r,
"z": z,
"c": c
}
#==================================================================
#
......@@ -84,14 +112,21 @@ ssp_data = {"r":r,"z":z,"c":c}
#
#==================================================================
bunits = '''W''' # Bottom absorption units: (dB/m)*kHZ
itype = '''R''' # Not used, but required...
xbty = array([x1,x2])/1000.0
ybty = array([-8.0,8.0])
zbty = array([[z1,z1],[z2,z2]])
p = array([1700.0,0.0,1.99,0.5])
bottom_data = {"itype":itype,"xbty":xbty,"ybty":ybty,"zbty":zbty,"p":p,"units":bunits}
bunits = '''W''' # Bottom absorption units: (dB/m)*kHZ
itype = '''R''' # Not used, but required...
xbty = array([x1, x2])/1000.0
ybty = array([-8.0, 8.0])
zbty = array([[z1, z1], [z2, z2]])
p = array([1700.0, 0.0, 1.99, 0.5])
bottom_data = {
"itype": itype,
"xbty": xbty,
"ybty": ybty,
"zbty": zbty,
"p": p,
"units": bunits
}
#==================================================================
#
......@@ -99,8 +134,9 @@ bottom_data = {"itype":itype,"xbty":xbty,"ybty":ybty,"zbty":zbty,"p":p,"units":b
#
#==================================================================
yarraykm = arange(0.035,5.000,0.005); nya = yarraykm.size
zarray = array([9.3])
yarraykm = arange(0.035, 5.000, 0.005)
nya = yarraykm.size
zarray = array([9.3])
options1 = '''SVWT''' # No ati file expected
options2 = '''A*''' # bty file expected
......@@ -112,12 +148,20 @@ nr = nya
rd = zarray
r = array([min(yarraykm), max(yarraykm)])
options = {"options1":options1,"options2":options2,"options3":options3,"options4":options4,\
"nrd":nrd,"nr":nr,"rd":rd,"r":r}
options = {
"options1": options1,
"options2": options2,
"options3": options3,
"options4": options4,
"nrd": nrd,
"nr": nr,
"rd": rd,
"r": r
}
print("Writing environmental file...")
wbellhop3denv('wedge',case_title,source_data,surface_data,ssp_data,bottom_data,options)
wbellhop3denv('wedge', case_title, source_data, surface_data, ssp_data, bottom_data, options)
print( "Running Bellhop3D..." )
......@@ -126,7 +170,7 @@ system("bellhop3d.exe wedge")
print( "Reading output data..." )
filename = 'wedge.shd'
pressure,geometry = readshd(filename,xs,ys)
pressure, geometry = readshd(filename, xs, ys)
p = squeeze( pressure, axis=(0,1) )
p = squeeze( p )
......
......@@ -3,9 +3,12 @@
# Bellhop: Block problem
# Mexilhoeira Grande, Dom Dez 25 22:07:45 WET 2016
# Written by Tordar
#
# Modified by David Van Komen (BYU)
#
#=======================================================================
# TODO: Fix these messy imports
from os import system
from numpy import *
from scipy.io import *
......@@ -17,17 +20,17 @@ from plotray import *
case_title = "Block problem"
freq = 100.0
Rmaxkm = 5.0; Rmax = Rmaxkm*1000.0
Dmax = 1000.0
cw = 1500.0 # sound speed in water
cb = 1700.0 # sound speed in lower halfspace
rhow = 1.0 # density in water
rhob = 2.0 # density in lower halfspace
freq = 100.0
Rmaxkm = 5.0; Rmax = Rmaxkm*1000.0 # converts from km to m
Dmax = 1000.0
cw = 1500.0 # sound speed in water
cb = 1700.0 # sound speed in lower halfspace
rhow = 1.0 # density in water
rhob = 2.0 # density in lower halfspace
source_nrays = 301 # number of propagation rays considered #
source_aperture = 20.0 # maximum launching angle (degrees) #
source_ray_step = 5.0 # ray calculation step (meters) #
source_nrays = 301 # number of propagation rays considered #
source_aperture = 20.0 # maximum launching angle (degrees) #
source_ray_step = 5.0 # ray calculation step (meters) #
#==================================================================
#
......@@ -35,17 +38,24 @@ source_ray_step = 5.0 # ray calculation step (meters) #
#
#==================================================================
nzs = 1
zs = array([500.0])
rs = array([ 0.0])
zbox = 1001.0
rbox = 5.1 # km!!!!!
box = array([source_ray_step,zbox,rbox])
thetas = array([source_nrays,-source_aperture,source_aperture])
p = zeros(1)
comp = ''
source_data = {"zs":zs, "box":box, "f":freq, "thetas":thetas, "p":p, "comp":comp}
nzs = 1
zs = array([500.0])
rs = array([ 0.0])
zbox = 1001.0
rbox = 5.1 # km!!!!!
box = array([source_ray_step,zbox,rbox])
thetas = array([source_nrays,-source_aperture,source_aperture])
p = zeros(1)
comp = ''
source_data = {
"zs": zs,
"box": box,
"f": freq,
"thetas": thetas,
"p": p,
"comp": comp
}
#==================================================================
#
......@@ -53,12 +63,17 @@ source_data = {"zs":zs, "box":box, "f":freq, "thetas":thetas, "p":p, "comp":comp
#
#==================================================================
itype = ''
xati = [] # The *.ati file won't be written
p = [] # Surface properties
aunits= ''
itype = ''
xati = [] # The *.ati file won't be written
p = [] # Surface properties
aunits = ''
surface_data = {"itype":itype,"x":xati,"p":p,"units":aunits}
surface_data = {
"itype": itype,
"x": xati,
"p": p,
"units": aunits
}
#==================================================================
#
......@@ -70,7 +85,11 @@ z = array([0.0,Dmax])
c = array([cw,cw])
r = []
ssp_data = {"r":r,"z":z,"c":c}
ssp_data = {
"r":r,
"z":z,
"c":c
}
#==================================================================
#
......@@ -78,14 +97,20 @@ ssp_data = {"r":r,"z":z,"c":c}
#
#==================================================================
rbty = array([rs[0]-2,2000,2010,2990,3000,Rmax+2]); rbtykm = rbty/1000.0
zbty = array([Dmax,Dmax,500,500,Dmax,Dmax])
itype = '''L''' # RID properties
rbty = array([rs[0]-2, 2000, 2010, 2990, 3000, Rmax+2])
rbtykm = rbty/1000.0 # km conversion
zbty = array([Dmax, Dmax, 500, 500, Dmax, Dmax])
itype = '''L''' # RID properties
bunits = '''W'''
xbty = array([rbtykm,zbty]) # Bottom coordinates
p = array([2000.0,0.0,2.0,0.5,0.0]) # Bottom properties
xbty = array([rbtykm, zbty]) # Bottom coordinates
p = array([2000.0, 0.0, 2.0, 0.5, 0.0]) # Bottom properties
bottom_data = {"itype":itype,"x":xbty,"p":p,"units":bunits}
bottom_data = {
"itype": itype,
"x": xbty,
"p": p,
"units": bunits
}
#==================================================================
#
......@@ -98,24 +123,35 @@ options2 = '''A*'''
options3 = '''R'''
options4 = []
rarray = zeros(1); rarraykm = zeros(1)
zarray = zeros(1)
rarray = zeros(1)
rarraykm = zeros(1)
zarray = zeros(1)
options = {"options1":options1,"options2":options2,"options3":options3,"options4":options4,"rarray":rarraykm,"zarray":zarray}
options = {
"options1": options1,
"options2": options2,
"options3": options3,
"options4": options4,
"rarray": rarraykm,
"zarray": zarray
}
print("Writing environmental file...")
wbellhopenvfil('block',case_title,source_data,surface_data,ssp_data,bottom_data,options)
wbellhopenvfil('block', case_title, source_data, surface_data, ssp_data,
bottom_data, options)
print( "Running Bellhop..." )
# TODO: WHAT??? It apparantly needs a bellhop.exe to work????
# That's something I might need to build and even test
system("bellhop.exe block")
print( "Reading output data..." )
figure(1)
plotray('block.ray')
plot(rbty,-zbty,'k',linewidth=2)
plot(rbty, -zbty, 'k', linewidth=2)
xlabel('Range (m)')
ylabel('Depth (m)')
title('Bellhop - Block problem')
......@@ -123,4 +159,3 @@ title('Bellhop - Block problem')
show()
print("done.")
......@@ -4,8 +4,12 @@
# Mexilhoeira Grande, Dom Dez 25 22:05:53 WET 2016
# Written by Tordar
#
# Modified by David Van Komen (BYU)
#
#=======================================================================
# TODO: Fix these messy imports
# That also means fixing all of the calls to functions...
from os import system
from numpy import *
from scipy.io import *
......@@ -17,17 +21,18 @@ from readshd import *
case_title = "Block problem"
freq = 100.0
Rmaxkm = 5.0; Rmax = Rmaxkm*1000.0
Dmax = 1000.0
cw = 1500.0 # sound speed in water
cb = 1700.0 # sound speed in lower halfspace
rhow = 1.0 # density in water
rhob = 2.0 # density in lower halfspace
freq = 100.0
Rmaxkm = 5.0
Rmax = Rmaxkm*1000.0
Dmax = 1000.0
cw = 1500.0 # sound speed in water
cb = 1700.0 # sound speed in lower halfspace
rhow = 1.0 # density in water
rhob = 2.0 # density in lower halfspace
source_nrays = 301 # number of propagation rays considered #
source_aperture = 20.0 # maximum launching angle (degrees) #
source_ray_step = 5.0 # ray calculation step (meters) #
source_nrays = 301 # number of propagation rays considered #
source_aperture = 20.0 # maximum launching angle (degrees) #
source_ray_step = 5.0 # ray calculation step (meters) #
#==================================================================
#
......@@ -35,17 +40,24 @@ source_ray_step = 5.0 # ray calculation step (meters) #
#
#==================================================================
nzs = 1
zs = array([500.0])
rs = array([ 0.0])
zbox = 1001.0
rbox = 5.1 # km!!!!!
box = array([source_ray_step,zbox,rbox])
thetas = array([source_nrays,-source_aperture,source_aperture])
p = zeros(1)
comp = ''
source_data = {"zs":zs, "box":box, "f":freq, "thetas":thetas, "p":p, "comp":comp}
nzs = 1
zs = array([500.0])
rs = array([ 0.0])
zbox = 1001.0
rbox = 5.1 # km!!!!!
box = array([source_ray_step, zbox, rbox])
thetas = array([source_nrays, -source_aperture, source_aperture])
p = zeros(1)
comp = ''
source_data = {
"zs": zs,
"box": box,
"f": freq,
"thetas": thetas,
"p": p,
"comp": comp
}
#==================================================================
#
......@@ -53,12 +65,17 @@ source_data = {"zs":zs, "box":box, "f":freq, "thetas":thetas, "p":p, "comp":comp
#
#==================================================================
itype = ''
xati = [] # The *.ati file won't be written
p = [] # Surface properties
aunits= ''
itype = ''
xati = [] # The *.ati file won't be written
p = [] # Surface properties
aunits = ''
surface_data = {"itype":itype,"x":xati,"p":p,"units":aunits}
surface_data = {
"itype": itype,
"x": xati,
"p": p,
"units": aunits
}
#==================================================================
#
......@@ -66,11 +83,15 @@ surface_data = {"itype":itype,"x":xati,"p":p,"units":aunits}
#
#==================================================================
z = array([0.0,Dmax])
c = array([cw,cw])
z = array([0.0, Dmax])
c = array([cw, cw])
r = []
ssp_data = {"r":r,"z":z,"c":c}
ssp_data = {
"r": r,
"z": z,
"c": c
}
#==================================================================
#
......@@ -78,14 +99,20 @@ ssp_data = {"r":r,"z":z,"c":c}
#
#==================================================================
rbty = array([rs[0]-2,2000,2010,2990,3000,Rmax+2]); rbtykm = rbty/1000.0
zbty = array([Dmax,Dmax,500,500,Dmax,Dmax])
rbty = array([rs[0]-2, 2000, 2010, 2990, 3000, Rmax+2])
rbtykm = rbty/1000.0
zbty = array([Dmax, Dmax, 500, 500,Dmax, Dmax])
itype = '''L''' # RID properties
bunits = '''W'''
xbty = array([rbtykm,zbty]) # Bottom coordinates
p = array([2000.0,0.0,2.0,0.5,0.0]) # Bottom properties
xbty = array([rbtykm, zbty]) # Bottom coordinates
p = array([2000.0, 0.0, 2.0, 0.5, 0.0]) # Bottom properties
bottom_data = {"itype":itype,"x":xbty,"p":p,"units":bunits}
bottom_data = {
"itype": itype,
"x": xbty,
"p": p,
"units":bunits
}
#==================================================================
#
......@@ -98,14 +125,22 @@ options2 = '''A*'''
options3 = '''C''' # Coherent TL (rewrite final block accordingly)
options4 = []
rarray = linspace(0,Rmax,501); rarraykm = rarray/1000.0
zarray = linspace(0,Dmax,101)
rarray = linspace(0, Rmax, 501)
rarraykm = rarray/1000.0
zarray = linspace(0, Dmax, 101)
options = {"options1":options1,"options2":options2,"options3":options3,"options4":options4,"rarray":rarraykm,"zarray":zarray}
options = {
"options1": options1,
"options2": options2,
"options3": options3,
"options4": options4,
"rarray": rarraykm,
"zarray": zarray
}
print("Writing environmental file...")
wbellhopenvfil('block',case_title,source_data,surface_data,ssp_data,bottom_data,options)
wbellhopenvfil('block', case_title, source_data, surface_data, ssp_data, bottom_data, options)