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.

It’s been a long time since I’ve

(a) used Python, numpy, scipy, or matplotlib

(b) done a beam pattern calculation of any kind

So I thought it would be a good chance to brush up on several things all at once, and maybe even have the satisfaction of making a pretty picture at the end of the night.

Well, here’s the (perhaps not so) pretty picture:

This was so simple, I’m almost embarrassed to post it, but it will probably be helpful for me later, so here it is.  What is it?  I pretended I had two little pingers spaced 5 wavelengths apart, sending out continuous sinusoidal waves.  I summed these signals at each of the elements in a 2-D grid, plotted the resulting magnitudes, and voila.

What I had to learn:

- How to use exp() using numpy – I had to load it separately for some reason, otherwise it would default to the exp() that comes with the basic math package, which doesn’t allow computations on arrays.

-How to plot an image, a la “imagesc” in Matlab.  This is accomplished using <code>imshow(datahere)</code>, and a colorbar can be added using <code>colorbar()</code>.

-How to do a screen capture of a certain area on a Mac computer:  command+shift+4 gives you crosshairs so you can pick what you like.

It would have been nice to add some attenuation, but I’m too tired now.  One more thing that I need to figure out though, what’s the best way to post my code?  I probably should figure out how to put it up on the bluehost server (where this blog lives), and link to it from here.

Wow, it was so weird to power up the new Macbook Pro.  It’s new, but still so familiar – like coming home after being away for a really long time.  It’s been nearly two years, after all, and ALL Linux during that time (at home, anyway).  Not to worry – I have not abandoned my Ubuntu computer!

So first things first – make sure I have what I need to dive back into Python.  Well, luckily, Python is built in already (2.6.1).  So there’s a step I don’t need to worry about.  Next up – Numpy and Scipy.  It’s been so long since I’ve done this.  When installing open source packages on my Mac before, I used Fink.  I have vague memories of Val and Kurt walking me through it way back when.  Unfortunately I didn’t document it very well the first time.  I’ll keep better track this time.

Some notes on Fink install:

1) No binary for OS 10.6.x.  So I need to do a source install.  Instructions here, starting with downloading the source tarball.

2)  Install Xcode Tools (required for installing packages from source).

• registered as an Apple Developer (woo, fancy), and downloaded xcode_3.2.3_and_ios_sdk_4.0.1.dmg (good lord, it’s 2.28GB!)
• Run xcode installation

3) unpack the fink tarball

4) Install Fink base setup using ./bootstrap

5) Set PATH environment variable using /sw/bin/pathsetup.sh

And there’s Fink, all set up (hopefully).  It’s getting late, I’ll save Numpy and Scipy for another day.

I just tried running David’s script, and a waterfall of numbers flew by the terminal window.  Too many – I couldn’t scroll up to the beginning to see if the script was doing what I thought it was.  So how to direct output to a file?  A quick Google search brought me to this site, and below is the useful part:

< Standard output >

Many Linux commands print their output to screen. For example, ls does this when it lists the contents of a directory: you see the output, the directory listing, on your screen.cat does the same: it concatenates a file and sends the results to your screen, and thus you can see the file’s contents. But the screen isn’t the only place where the commands can print their output because you can redirect the output of several commands to files, devices, and even to the input of other commands.

The CLI programs that display their results do so usually by sending the results to standard output, or stdout for short. By default, standard output directs its contents to the screen, as you’ve seen with the ls and cat commands. But if you want to direct the output to somewhere else, you can use the > character. For example, to redirect the output to a file, you can use the > character like this:
$ ls > dir_listing.txt

The above redirects the output of the ls command to a file called dir_listing.txt. Because the output is redirected to a file, you don’t see any results of ls on your screen.

Each time you repeat the above command, the file dir_listing.txt is overwritten with the results of the ls command. If you want to append the new results to the file instead of rewriting it, you can use >> instead:
$ ls >> dir_listing.txt

Each time you repeat the above command, the new output of ls is added at the end of the dir_listing.txt file instead of overwriting the file.

The following adds the contents of File1 at the end of File2:
$ cat File1 >> File2

Like I told you before, files aren’t the only places where you can redirect the standard output. You can redirect it to devices, too:
$ cat sound.wav > /dev/audio

As you saw, in the above example the cat command concatenates a file named sound.wav and the results are sent to a device called /dev/audio. What’s the fun here, then?/dev/audio is the audio device, and when you send there the contents of a sound file, that file is played. So if your sound is properly configured, the above command plays the file sound.wav!

Once I did that, I was able to save the file, and I could see that it was, indeed, doing what it was supposed to be doing (of course it was – it was only a matter of me figuring out how to look at it!).

Here’s the top part of the output (before it turns into a waterfall of numbers):

Bytes: 396 RecordType: 7200 fragmentNumber: 0
Bytes: 4613 RecordType: 7001 fragmentNumber: 0
Bytes: 4837 RecordType: 7000 fragmentNumber: 0
Bytes: 13109 RecordType: 7004 fragmentNumber: 0
Bytes: 21905 RecordType: 7006 fragmentNumber: 0
Bytes: 33072 RecordType: 7027 fragmentNumber: 0
%% DATA RECORD FRAME (11167) bytes
ProtocalVersion: 5
Offset: 60
Sync Pattern: 65535
Size: 11167
Optional Data Offset: 0
Optional Data Identifier: 0
Year: 2009
Day: 353
Second: 55.214169
Hour: 17
Minute: 53
Record Type Identifier: 7027
Device Identifier: 7125
System Enumerator: 0
Flags: 32769
Total Records In Fragmented Data Record Set: 0
Fragment Number: 0
DATA SECTION: 11099 (bytes)
Checksum: 788564
--- Data Section --
(Detection Properties Record)
Sonar Id: 0
Ping Number: 0
Multi-Ping Sequence: 0
Points(?): 500
Size(?): 22
Sample Rate: 34482.757812
Transducer Angle: 0.000000

It’s sonar data! It really is! Okay, mostly header info, but still… It’s probably a bit silly to be this excited about simply running someone else’s program. But I had to figure out at least 4 or 5 things just to get to that point, so it’s not a wasted evening, right?  :-)

Hey – it’s been a while…but I’m back from Peru, ready to jump back on the Linux/Ubuntu/Python/GMT/MBsystem train – especially since I’ve gotten some tips in the meantime that will hopefully get me back on track with the issues I was having.  (I’m talking about you, GDAL! Thanks Kurt and Monica for your help!)

Happy Pi day, everyone!