Today was my last day at Reson, and I’m finally getting around to writing the Python version of my ray tracing script. What I want to do is take sound speed profiles as input, and compute the difference in depth that would result from using one versus the other at different multibeam launch angles. It’s just a little tool that I’d like to have to figure out how much error I could expect given two profiles.
So far, I’ve read in the profile (from a text file), and parsed the data:
filename = 'ssp.txt'
ssp = []
depth = []
file1 = open(filename)
for line in file1.readlines():
# print line
splitup = line.split()
ssp.append(float(splitup[1]))
depth.append(float(splitup[0]))
After this, I defined an array of launch angles (relative to horizontal, so +/- 75 degrees was actually 15 – 165), and an array of sound speed values interpolated from the sound speed profile to a certain depth increment (1m), and down to a given water depth (10m).
LaunchAngle = r_[15.:166.:5.] # array of launch angles in degrees
LaunchAngleRadians = LaunchAngle*pi/180
nadirdepth = 100 # nadir depth in meters
DepthIncrement = 0.5 # depth increment in m
if nadirdepth > max(depth):
ssp.append(0.016*(nadirdepth - max(depth)))
depth.append(nadirdepth)
if min(depth)>0:
ssp = insert(ssp,0,ssp[0])
depth = insert(depth,0,0)
depthvector = r_[0:nadirdepth:DepthIncrement]
interpfun = interp1d(depth,ssp)
sspvector = interpfun(depthvector)
And as if that much wasn’t suspicious enough, I decided that I wanted to do all of the calcs using matrix math. I figured I’d ditch all the loops that were in my original Matlab script, and replace them with some nice streamlined matrices. So efficient. The problem is that there are some exceptions that I need to take care of, that are difficult to deal with when doing everything in a series of matrices. For example, when the sound speed in a certain layer is the same as the sound speed in the first layer – then the gradient is zero, and the radius of curvature is infinite. Here’s what I’ve hacked together so far. It’s not pretty.
thetaLaunch = tile(LaunchAngleRadians,(size(depthvector),1))
depthmat = transpose(tile(depthvector,(size(LaunchAngleRadians),1)))
sspmat = transpose(tile(sspvector,(size(LaunchAngleRadians),1)))
g = (sspmat - sspvector[0]) / DepthIncrement
Rc = -sspvector[0]/(g*cos(thetaLaunch))
thetamat = arccos( ((Rc*cos(thetaLaunch)) - DepthIncrement ) / Rc )
dx = (Rc*sin(thetamat) - Rc*sin(thetaLaunch))
dxnandex = isnan(dx)
dx[dxnandex] = DepthIncrement/tan(thetaLaunch[dxnandex])
dx[0,] = 0
dxsum = cumsum(dx,axis=0)
Yuck. But at least I managed to get some rays that made sense.
The next step is to convert the ray paths to two-way travel times, and then compare the depth that would be calculated if the ray tracing were done with a second sound speed profile.
