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.

Thanks for the tip, Kurt!!

I guess this is probably a painfully obvious solution to anyone who knows what they’re doing, but it turns out I was missing the link to GDAL in my install_makefiles script. There’s a line in there that looks like this:


$LFLAGS = "-L$NETCDFLIBDIR -lm -lnetcdf";


It should have had a link to GDAL like this:


$LFLAGS = "-L$NETCDFLIBDIR -lm -lnetcdf -lgdal";


So it worked! It ran with no errors, and seemed to build correctly. So one step closer! I then continued along the instructions to setting the .bashrc file. I had left out one line before, so I added it, and then ran:


$ source .bashrc


To test my GMT and MBSystem installs, I tried:


$ psxy
$ man psxy
$ mbgrid
$ man mbgrid


Since I don’t know how to use GMT and MBSystem yet, I don’t know if these did what they were supposed to. But it seemed like they did, except for mbgrid, which came up with an error: “Unable to open data list file: datalist.mb-1“. But at least it recognized mbgrid as a command and tried to do something!

Ah, finally, I’m getting somewhere!

Thanks, Kurt, for the tip.  Here’s a link to Larry’s interview with CNN, where he talks the search for “phantom” oil plumes.

CNN Video – Phantom plumes

Both of my master’s thesis advisors were down in the Gulf of Mexico recently, searching for these plumes using various types of sensors.

Kurt also linked to the recently released cruise report for the mid-water mapping mission in the Gulf of Mexico:  NOAA Ship Thomas Jefferson Deepwater Horizon Response Mission Report

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?  :-)

Jonathan Beaudoin, an old friend of mine from UNB, has recently moved down to UNH… and he’s got a blog!  Yay for nerdy blogs!  Check it out here:  2bitbrain.blogspot.com.  It’s totally on my Google Reader list, as of 5 minutes ago.

Jonathan does some really neat research, and is a fantastic combination of computer scientist/hydrographer /general smarty-pants, so I’m looking forward to following his adventures as he moves into mac-land.