I like Sonars. I like computers. I like science. I'm a geek

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.

I got some Fink tips from Kurt on my post from a couple of days ago. I’m going to re-post them here:

fink configure # Say yes to unstable (which means gets actual updates – stable is missing lots)
fink selfupdate-rsync (do this once to get rsync, after just do fink selfupdate)
fink list py26 # See all the fink pythkn 2.6 packages

Once I’d done this, I installed Scipy using:

fink install scipy-py26

Numpy was installed automatically with scipy.

Next up, matplotlib:

fink install matplotlib-py26

As part of a little script I’m writing, I need to do some simple linear interpolation. The Matlab equivalent of what I’m doing is called ‘interp1′. So far I’ve come across two ways to do 1-d, linear interpolation. One is ‘interp’ in numpy. The other is ‘interp1d’ in scipy.interpolate. I’m not completely sure of the difference. Is one faster? Is one older? I’ve googled around a bit, but still haven’t found a clear answer. For now, I’ve implemented ‘interp1d’, which is less similar to Matlab’s ‘interp1′. You first define an interpolated object given your x and y vectors. Then you call that object with your new x-values to generate the new y-values.

And here is an example, chopped out of my code:

depthvector1 = r_[0:nadirdepth:DepthIncrement]
interpfun1 = interp1d(depth1,ssp1) # my x and y vectors
sspvector1 = interpfun1(depthvector1) # find the new y's for this x

Sometimes I think I have my biggest breakthroughs while biking home. Tonight, for instance, it dawned on me (on the stretch between the beach and Patterson) just how awful the Matlab version of my ray tracing code really is. And I realised that I could vastly improve its simplicity and efficiency by using matrices (Yikes! Why on earth did I implement all of those nested loops?). I pictured the steps and the math involved, the Matlab commands that I would need, and got really excited – until I realized that I don’t know how to do any of the same stuff in Python yet. First of all, I need to figure out how to work with arrays in Matlab. Thankfully, there’s a special page on Scipy.org Matlab users who are trying to migrate to Numpy/Scipy. Yes! And here it is: Numpy for Matlab Users. My first question was, “How do I make it do element-wise operations?” In Matlab, you need to modify the operator, for example using “.*” to multiply by elements, and just “*” to do matrix multiplication. Matlab thinks in linear algebra for 2d matrices. But, as I quickly learned, the default array in Numpy does element by element operations, unless you create a special matrix object. Awesome.

I just picked up my second python book – although this is probably the one I should have bought first.  It’s ‘Learning Python’ from O’Reilly.  I already had Beginning Python.  I’ve been wanting to get the O’Reilly one for a while now, but never really wanted to spend the 50 bucks.  I finally broke down and got it at Borders this morning.  I hope it’s a good reference.  The first thing I’d like to do is to remember how to do things like read files, parse data, create plots.  I think I remember having to add Scipy and Matplotlib.

If I want Scipy, I need to also get Numpy.  In the Scipy/Numpy installation instructions for Ubuntu, it gives a couple of options for installation.  One of these is from a guy called Andrew Straw, who has an unofficial repository for NumPy .deb packages.  I’m going to try that first.  I just checked out his repository and he just updated for Karmic on 2 December.  I’m using the instructions here.