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.

Here we go again… I was trying to draw a little diagram, and since I’m using a Mac, figured I’d look up OmniGraffle.  But OmniGraffle is now $100 for the basic package – and $200 for professional.  It wouldn’t be quite so painful except I got it for free back when I got my old Mac Powerbook back in 2006.

So it’s back to Inkscape.  I love open source!  It was an easy install – I decided to take a pass on compiling from source, and just grabbed the .dmg file from the downloads page. And of course, I already went through the whole Inkscape plus Latex thing a while ago.

And, for fun, here’s a screen shot of Inkscape in action:

… who needs OmniGraffle anyway… (sniffle)

Yeah – thanks, Val!

Square Dancing Article from The New York Times

This article made me smile – it gives a fun description of the Pythagorean theorem. Yay, geometry! My favorite is the part where the author tells where the word geometry comes from geo (earth) and metry (measurement) – geometry was originally used in solving problems of land measurement. This makes me happy to be a surveyor

Okay, I’m being really indecisive here… but I finally found one that I think I like (so far).  It’s not perfect, and it’s a bit quirky, but at least I managed to do what I need to do.  Here’s the main website:  http://www.elisanet.fi/ptvirtan/software/textext/. I also found this blog post to be helpful.

I tried making a really simple depiction of Snell’s law, basically showing the relationship between the angles of incidence and refraction and the refractive indices of two layers. In this case, it is sound waves travelling through homogeneous layers of water with different sound speeds. Easy peasy The sketch only took me about 15 minutes to create.

All I had to do was go to Extensions –> textext, set the scale to 2.5, and enter the following text into the box:

\begin{minipage}{5cm}

$\frac{c_{i}}{\sin{\theta_{i}}} = \frac{c_{i+1}}{\sin{\theta_{i+1}}}$

\end{minipage}

Yes, it’s sort of annoying to have to add that \begin{minipage} stuff. I’d much rather have it more behind the scenes, but I suppose that does give you more control. The resulting image looked like this:

So, nothing too special, but it’ll work for now.  I’ll keep an eye out for other options, but I am getting to really like Inkscape, so maybe I’ll stick with it.  There might be other Latex plugins too that work better.  This one is convenient because you can re-edit an equation once it’s been entered.  If it’s a complicated equation, it can be a real pain to re-type the whole thing if there’s only one little mistake.

One more thing I ought to mention: I keep getting this error message. I’ve been basically ignoring it because I’m getting the correct output, but Inkscape is clearly not happy with something:

/home/michelle/.config/inkscape/extensions/textext.py:55: DeprecationWarning: the md5 module is deprecated; use hashlib instead
import os, sys, tempfile, traceback, glob, re, md5, copy

PS: Why did I abandon SketchLatex so quickly, when it seemed so promising this morning? Simply because their website was down. Anyway, the Skencil GUI was a bit old-fashioned Unix-y anyway. That alone would not have stopped me from using it, but combined with other difficulties it becomes a problem. Hence going back to Inkscape
PPS: The only Inkscape drawing I’ve ever done was my little star website icon (see it up there? the purple star on the blue background?)

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.