-
Notifications
You must be signed in to change notification settings - Fork 1
/
bellhop_block_tl.py
135 lines (105 loc) · 3.54 KB
/
bellhop_block_tl.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#=======================================================================
#
# Bellhop: Block problem
# Gambelas, qui 20 jun 2024 11:46:02
# Written by Orlando Camargo Rodriguez
#
#=======================================================================
from os import system
from numpy import *
from matplotlib.pyplot import *
from wbellhopenvfil import *
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
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 properties
#
#==================================================================
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}
#==================================================================
#
# Surface definition:
#
#==================================================================
itype = ''
xati = [] # The *.ati file won't be written
p = [] # Surface properties
aunits= ''
surface_data = {"itype":itype,"x":xati,"p":p,"units":aunits}
#==================================================================
#
# Sound speed:
#
#==================================================================
z = array([0.0,Dmax])
c = array([cw,cw])
r = []
ssp_data = {"r":r,"z":z,"c":c}
#==================================================================
#
# Bottom:
#
#==================================================================
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
bottom_data = {"itype":itype,"x":xbty,"p":p,"units":bunits}
#==================================================================
#
# Array:
#
#==================================================================
options1 = '''CVW''' # No ati file expected
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)
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)
print( "Running Bellhop..." )
system("bellhop.exe block")
print( "Reading output data..." )
filename = 'block.shd'
xs = nan
ys = nan
pressure,geometry = readshd(filename,xs,ys,freq)
p = squeeze( pressure, axis=(0,1) )
p = where( p == 0, nan, p )
tl = -20*log10( abs( p ) )
figure(1)
pcolormesh(rarray,zarray,tl,vmin=60,vmax=90,cmap='jet_r',shading='auto')
fill_between(rbty,zbty,Dmax)
colorbar()
plot(rs[0],zs[0],marker="<",markersize=16,color="k")
xlim(0,Rmax)
ylim(Dmax,0)
xlabel('Range (m)')
ylabel('Depth (m)')
title('Bellhop - Block problem')
show()
print("done.")