org-arl
Bellhop Arrivals Complex amplitude calculation
I believe I discovered an issue with arlpy's arrival amplitude calculation. When I change the frequency, the reported amplitude does not change for the identical beams/rays. After confirming the issue was not with the bellhop call, the env file, or the resulting arrival file, and in fact persists running arlpy functions against .arr files I produce using the test cases that came with bellhop, I found that this is due to not accounting for the complex part of the arrival delay when reporting arrival amplitude and time of arrival.
Lines 778-793 in uwapm.py currently reads:
data = self._readf(f, (float, float, float, float, float, float, int, int)) arrivals.append(_pd.DataFrame({ 'tx_depth_ndx': [j], 'rx_depth_ndx': [k], 'rx_range_ndx': [m], 'tx_depth': [tx_depth[j]], 'rx_depth': [rx_depth[k]], 'rx_range': [rx_range[m]], 'arrival_number': [n], 'arrival_amplitude': [data[0]_np.exp(1jdata[1])], 'time_of_arrival': [data[2]], 'angle_of_departure': [data[4]], 'angle_of_arrival': [data[5]], 'surface_bounces': [data[6]], 'bottom_bounces': [data[7]] }, index=[len(arrivals)+1]))
This does not account for field data[3], which is the complex part of the delay. When using the flag "T" in OPTION1 of the env file, this number is non-zero and accounts for the attenuation difference with frequency. If I look at the matlab code, read_arrivals_asc.m, I see that this field is handled by reporting the delay as a complex number: Arr.A( ir, 1:Narr, ird, isd ) = da( 1, 1:Narr ) .* exp( 1i * da( 2, 1:Narr ) * pi/180); Arr.delay( ir, 1:Narr, ird, isd ) = da( 3, 1:Narr ) + 1i * da( 4, 1:Narr ); Arr.SrcAngle( ir, 1:Narr, ird, isd ) = da( 5, 1:Narr ); Arr.RcvrAngle( ir, 1:Narr, ird, isd ) = da( 6, 1:Narr ); Arr.NumTopBnc( ir, 1:Narr, ird, isd ) = da( 7, 1:Narr ); Arr.NumBotBnc( ir, 1:Narr, ird, isd ) = da( 8, 1:Narr );
To fix the issue with arlpy.uwapm for arrivals files, there are two choices: 1) report delay/TOA as a complex number or 2) account for the complex delay and arrival phase shift in calculating of the arrival amplitude. I suggest the second, as it means you don't change the type entering the pandas data structure. In the below solution, I also fix the issue that I beleive is occuring due to not multiplying data[1] by pi/180 for the phase; that part of the arrivals file is reported in degrees, substantiated by the fact that it is always 180 degrees for a vacuum surface bounce and 0 for a direct path.
Here is what I did in my local copy (again Lines 778-793 in uwapm.py):
data = self._readf(f, (float, float, float, float, float, float, int, int)) arrivals.append(_pd.DataFrame({ 'tx_depth_ndx': [j], 'rx_depth_ndx': [k], 'rx_range_ndx': [m], 'tx_depth': [tx_depth[j]], 'rx_depth': [rx_depth[k]], 'rx_range': [rx_range[m]], 'arrival_number': [n], 'arrival_amplitude': [data[0]_np.exp(1j(data[1]math.pi/180+freqmath.pidata[3]1j+freqmath.pidata[2]))], 'time_of_arrival': [data[2]], 'angle_of_departure': [data[4]], 'angle_of_arrival': [data[5]], 'surface_bounces': [data[6]], 'bottom_bounces': [data[7]] }, index=[len(arrivals)+1])) :
This creates an arrival amplitude that accounts for phase shift due to propogation through the medium as well as the amplitude change indicated in the complex part of the delay. For now I am using this code in my local copy, but thought I should report the issue so that it can be resolved for others.
Please let me know if I can be helpful in understanding the issue, and thank you for providing/supporting this great acoustics modelling and signal processing resource!
Thanks, Erin Fischell, PhD Senior Scientist JPAnalytics, LLC
Thanks Erin. I agree with you that the second approach would be a better one. I'm fine with the change. Would you prefer to do a PR to get this in, since you've already done the work to implement?
@efischell Just wanted to touch base to know if you'll be raising a PR or should I be manually putting the changes in?
Since there is no PR on this, I can try to manually update and test as suggested above. These couple of weeks are busy, and so I'll only be able to do it end of the month. Alternatively, @tfabbri, if you want to take a stab at it, I'd be happy to review a PR.
@mchitre I will try to look at this issue and push a new PR. I'll keep you posted on it.
I have a worked a little on the report from @efischell.
Looking at the matlab code from read_arrivals_asc.m (latest version from April 2022) we have:
Arr( irr, irz, isd ).A = single(da( 1, 1 : Narr ) .* exp( 1.0i * da( 2, 1 : Narr ) * pi / 180.0 ) );
Arr( irr, irz, isd ).delay = single(da( 3, 1 : Narr ) + 1.0i * da( 4, 1 : Narr ) );
compared with the arlpy current reader a direct translation of the method would be:
arrivals.append(_pd.DataFrame({
...
'arrival_amplitude': [data[0]*_np.exp(1j*_np.deg2rad(data[1]))],
'time_of_arrival': [data[2] + 1j*data[3]],
...
}, index=[len(arrivals)+1]))
Considering the complex delay and arrival phase shift in calculating of the arrival amplitude as suggested @efischell we have a form $A e^{j(\phi + \omega t)}$ where $\phi=$data[1] and $\omega t =$freq * 2 * _np.pi *(data[2] + 1j* data[3]).
'arrival_amplitude': [data[0] * _np.exp( 1j * (_np.deg2rad(data[1]) + freq * 2 * _np.pi * (data[2] + 1j* data[3])))]
'time_of_arrival': [data[2] + 1j*data[3]],
Waiting your feedback before going through the pull-request.
Sounds good @tfabbri, for the arrival_amplitude. But the time_of_arrval should stay as data[2], since the complex part is accounted for in the amplitude change?
@mchitre , I agree on your point. The drawback is that we lose track of this data (complex part of time_of_arrival) if required by the user. I do not know which the best choice is.
For completeness, I would add some documentation to specify the difference between the arlpy reader and the matlab one.
Good idea.
So I can introduce with the next pull request a new field complex_time_of_arrival to keep track of this bellhop output.
'arrival_amplitude': [data[0] * _np.exp( 1j * (_np.deg2rad(data[1]) + freq * 2 * _np.pi * (data[2] + 1j* data[3])))]
'time_of_arrival': [data[2]],
'complex_time_of_arrival': [data[2] + 1j*data[3]],
'arrival_amplitude': [data[0] * _np.exp( 1j * (_np.deg2rad(data[1]) + freq * 2 * _np.pi * (data[2] + 1j* data[3])))]
this is creating issues with the newer AT models (released on 5/2023). As the above changes returns freq in tupple (please refer Traceback for more details).
can we do this :
# Extract the float value from the tuple
freq_value = freq[0]
# Use freq_value in your calculations
'arrival_amplitude': [data[0] * _np.exp(-1j * (_np.deg2rad(data[1]) + freq_value * 2 * _np.pi * (data[3] * 1j + data[2])))],
Here is the more detailed traceback :
2023-10-18 18:44:44.694 | ERROR | arlpy.uwapm:compute_arrivals:326 - An error has been caught in function 'compute_arrivals', process 'MainProcess' (32086), thread 'MainThread' (139911861442368):
Traceback (most recent call last):
File "/usr/lib/python3.8/runpy.py", line 194, in _run_module_as_main
return _run_code(code, main_globals, None,
│ │ └ {'__name__': '__main__', '__doc__': 'Entry point for launching an IPython kernel.\n\nThis is separate from the ipykernel pack...
│ └ <code object <module> at 0x7f3fc46200e0, file "/home/batman/.local/lib/python3.8/site-packages/ipykernel_launcher.py", line 1>
└ <function _run_code at 0x7f3fc45940d0>
File "/usr/lib/python3.8/runpy.py", line 87, in _run_code
exec(code, run_globals)
│ └ {'__name__': '__main__', '__doc__': 'Entry point for launching an IPython kernel.\n\nThis is separate from the ipykernel pack...
└ <code object <module> at 0x7f3fc46200e0, file "/home/batman/.local/lib/python3.8/site-packages/ipykernel_launcher.py", line 1>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel_launcher.py", line 17, in <module>
app.launch_new_instance()
│ └ <bound method Application.launch_instance of <class 'ipykernel.kernelapp.IPKernelApp'>>
└ <module 'ipykernel.kernelapp' from '/home/batman/.local/lib/python3.8/site-packages/ipykernel/kernelapp.py'>
File "/home/batman/.local/lib/python3.8/site-packages/traitlets/config/application.py", line 1053, in launch_instance
app.start()
│ └ <function IPKernelApp.start at 0x7f3fc068d310>
└ <ipykernel.kernelapp.IPKernelApp object at 0x7f3fc47174c0>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel/kernelapp.py", line 736, in start
self.io_loop.start()
│ │ └ <function BaseAsyncIOLoop.start at 0x7f3fc065da60>
│ └ <tornado.platform.asyncio.AsyncIOMainLoop object at 0x7f3fbe5f6700>
└ <ipykernel.kernelapp.IPKernelApp object at 0x7f3fc47174c0>
File "/home/batman/.local/lib/python3.8/site-packages/tornado/platform/asyncio.py", line 195, in start
self.asyncio_loop.run_forever()
│ │ └ <function BaseEventLoop.run_forever at 0x7f3fc23038b0>
│ └ <_UnixSelectorEventLoop running=True closed=False debug=False>
└ <tornado.platform.asyncio.AsyncIOMainLoop object at 0x7f3fbe5f6700>
File "/usr/lib/python3.8/asyncio/base_events.py", line 570, in run_forever
self._run_once()
│ └ <function BaseEventLoop._run_once at 0x7f3fc2306430>
└ <_UnixSelectorEventLoop running=True closed=False debug=False>
File "/usr/lib/python3.8/asyncio/base_events.py", line 1859, in _run_once
handle._run()
│ └ <function Handle._run at 0x7f3fc23321f0>
└ <Handle <TaskWakeupMethWrapper object at 0x7f3f91a0d250>(<Future finis...0f0>, ...],))>)>
File "/usr/lib/python3.8/asyncio/events.py", line 81, in _run
self._context.run(self._callback, *self._args)
│ │ │ │ │ └ <member '_args' of 'Handle' objects>
│ │ │ │ └ <Handle <TaskWakeupMethWrapper object at 0x7f3f91a0d250>(<Future finis...0f0>, ...],))>)>
│ │ │ └ <member '_callback' of 'Handle' objects>
│ │ └ <Handle <TaskWakeupMethWrapper object at 0x7f3f91a0d250>(<Future finis...0f0>, ...],))>)>
│ └ <member '_context' of 'Handle' objects>
└ <Handle <TaskWakeupMethWrapper object at 0x7f3f91a0d250>(<Future finis...0f0>, ...],))>)>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel/kernelbase.py", line 516, in dispatch_queue
await self.process_one()
│ └ <function Kernel.process_one at 0x7f3fc0baa700>
└ <ipykernel.ipkernel.IPythonKernel object at 0x7f3fbe5f6d60>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel/kernelbase.py", line 505, in process_one
await dispatch(*args)
│ └ ([<zmq.sugar.frame.Frame object at 0x7f3fbc5c1e00>, <zmq.sugar.frame.Frame object at 0x7f3fbc5c15c0>, <zmq.sugar.frame.Frame ...
└ <bound method Kernel.dispatch_shell of <ipykernel.ipkernel.IPythonKernel object at 0x7f3fbe5f6d60>>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel/kernelbase.py", line 412, in dispatch_shell
await result
└ <coroutine object Kernel.execute_request at 0x7f3fbc5dc340>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel/kernelbase.py", line 740, in execute_request
reply_content = await reply_content
└ <coroutine object IPythonKernel.do_execute at 0x7f3f91a18b40>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel/ipkernel.py", line 422, in do_execute
res = shell.run_cell(
│ └ <function ZMQInteractiveShell.run_cell at 0x7f3fc067baf0>
└ <ipykernel.zmqshell.ZMQInteractiveShell object at 0x7f3fbe602220>
File "/home/batman/.local/lib/python3.8/site-packages/ipykernel/zmqshell.py", line 546, in run_cell
return super().run_cell(*args, **kwargs)
│ └ {'store_history': True, 'silent': False, 'cell_id': 'eb28d754-2fbd-4ba5-bfe9-a44146ce96bb'}
└ ('arrivals = pm.compute_arrivals(env)\npm.plot_arrivals(arrivals, width=900)',)
File "/home/batman/.local/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3009, in run_cell
result = self._run_cell(
│ └ <function InteractiveShell._run_cell at 0x7f3fc15075e0>
└ <ipykernel.zmqshell.ZMQInteractiveShell object at 0x7f3fbe602220>
File "/home/batman/.local/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3064, in _run_cell
result = runner(coro)
│ └ <coroutine object InteractiveShell.run_cell_async at 0x7f3fbc5dc540>
└ <function _pseudo_sync_runner at 0x7f3fc14f0e50>
File "/home/batman/.local/lib/python3.8/site-packages/IPython/core/async_helpers.py", line 129, in _pseudo_sync_runner
coro.send(None)
│ └ <method 'send' of 'coroutine' objects>
└ <coroutine object InteractiveShell.run_cell_async at 0x7f3fbc5dc540>
File "/home/batman/.local/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3269, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
│ │ │ │ └ '/tmp/ipykernel_32086/1324318870.py'
│ │ │ └ [<_ast.Assign object at 0x7f3f8e15eb50>, <_ast.Expr object at 0x7f3f8e15ed90>]
│ │ └ <_ast.Module object at 0x7f3f8dd83040>
│ └ <function InteractiveShell.run_ast_nodes at 0x7f3fc15078b0>
└ <ipykernel.zmqshell.ZMQInteractiveShell object at 0x7f3fbe602220>
File "/home/batman/.local/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3448, in run_ast_nodes
if await self.run_code(code, result, async_=asy):
│ │ │ │ └ False
│ │ │ └ <ExecutionResult object at 7f3f8e15ec40, execution_count=5 error_before_exec=None error_in_exec=None info=<ExecutionInfo obje...
│ │ └ <code object <module> at 0x7f3f8e08c9d0, file "/tmp/ipykernel_32086/1324318870.py", line 1>
│ └ <function InteractiveShell.run_code at 0x7f3fc1507940>
└ <ipykernel.zmqshell.ZMQInteractiveShell object at 0x7f3fbe602220>
File "/home/batman/.local/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3508, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
│ │ │ │ └ {'__name__': '__main__', '__doc__': 'Automatically created module for IPython interactive environment', '__package__': None, ...
│ │ │ └ <ipykernel.zmqshell.ZMQInteractiveShell object at 0x7f3fbe602220>
│ │ └ <property object at 0x7f3fc14f75e0>
│ └ <ipykernel.zmqshell.ZMQInteractiveShell object at 0x7f3fbe602220>
└ <code object <module> at 0x7f3f8e08c9d0, file "/tmp/ipykernel_32086/1324318870.py", line 1>
File "/tmp/ipykernel_32086/1324318870.py", line 1, in <module>
arrivals = pm.compute_arrivals(env)
│ │ └ {'name': 'arlpy', 'type': '2D', 'frequency': 25000, 'soundspeed': 1500, 'soundspeed_interp': 'spline', 'bottom_soundspeed': 1...
│ └ <function compute_arrivals at 0x7f3f8e160310>
└ <module 'arlpy.uwapm' from '/home/batman/.local/lib/python3.8/site-packages/arlpy/uwapm.py'>
> File "/home/batman/.local/lib/python3.8/site-packages/arlpy/uwapm.py", line 326, in compute_arrivals
return model.run(env, arrivals, debug)
│ │ │ │ └ False
│ │ │ └ 'arrivals'
│ │ └ {'name': 'arlpy', 'type': '2D', 'frequency': 25000, 'soundspeed': 1500, 'soundspeed_interp': 'spline', 'bottom_soundspeed': 1...
│ └ <function _Bellhop.run at 0x7f3f8e160af0>
└ <arlpy.uwapm._Bellhop object at 0x7f3fbc5bbc70>
File "/home/batman/.local/lib/python3.8/site-packages/arlpy/uwapm.py", line 611, in run
results = taskmap[task][1](fname_base)
│ │ └ '/tmp/tmpa3krxug6'
│ └ 'arrivals'
└ {'arrivals': ['A', <bound method _Bellhop._load_arrivals of <arlpy.uwapm._Bellhop object at 0x7f3fbc5bbc70>>], 'eigenrays': [...
File "/home/batman/.local/lib/python3.8/site-packages/arlpy/uwapm.py", line 800, in _load_arrivals
'arrival_amplitude': [data[0] * _np.exp( -1j * (_np.deg2rad(data[1]) + freq * 2 * _np.pi * (data[3] * 1j + data[2])))],
│ │ │ │ │ │ │ │ │ │ └ (1.25387023e-05, 1636.07947, 0.721795559, -4.85707142e-06, -22.5382519, 22.5382519, 9, 8)
│ │ │ │ │ │ │ │ │ └ (1.25387023e-05, 1636.07947, 0.721795559, -4.85707142e-06, -22.5382519, 22.5382519, 9, 8)
│ │ │ │ │ │ │ │ └ 3.141592653589793
│ │ │ │ │ │ │ └ <module 'numpy' from '/home/batman/.local/lib/python3.8/site-packages/numpy/__init__.py'>
│ │ │ │ │ │ └ (25000.0,)
│ │ │ │ │ └ (1.25387023e-05, 1636.07947, 0.721795559, -4.85707142e-06, -22.5382519, 22.5382519, 9, 8)
│ │ │ │ └ <ufunc 'deg2rad'>
│ │ │ └ <module 'numpy' from '/home/batman/.local/lib/python3.8/site-packages/numpy/__init__.py'>
│ │ └ <ufunc 'exp'>
│ └ <module 'numpy' from '/home/batman/.local/lib/python3.8/site-packages/numpy/__init__.py'>
└ (1.25387023e-05, 1636.07947, 0.721795559, -4.85707142e-06, -22.5382519, 22.5382519, 9, 8)
TypeError: can't multiply sequence by non-int of type 'float'
This is the issue current code raises :
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[5], line 2
1 arrivals = pm.compute_arrivals(env)
----> 2 pm.plot_arrivals(arrivals, width=900)
File ~/.local/lib/python3.8/site-packages/arlpy/uwapm.py:444, in plot_arrivals(arrivals, dB, color, **kwargs)
430 def plot_arrivals(arrivals, dB=False, color='blue', **kwargs):
431 """Plots the arrival times and amplitudes.
432
433 :param arrivals: arrivals times (s) and coefficients
(...)
442 >>> pm.plot_arrivals(arrivals)
443 """
--> 444 t0 = min(arrivals.time_of_arrival)
445 t1 = max(arrivals.time_of_arrival)
446 oh = _plt.hold()
AttributeError: 'NoneType' object has no attribute 'time_of_arrival'
Thanks @patel999jay, open to a PR on this. Could you raise one with the suggested changes?
@mchitre : do you want me to raise PR with this solution or have something else in mind ?
Yes, if that solution fixes it.
I don't use Python for my analysis anymore (having moved to Julia), so this gets no testing from me nowadays. If you're using it and have a solution, I'd appreciate a PR.