The Santa Barbara Channel Experiment (SBCX)

Propagation Codes

Most acoustic propagation codes will compile and run unmodified under Linux.
Use f2c rather than g77 or egcs for compilation. g77 is still a beta compiler, and there are stability issues associated with it. f2c, when combined with the fort77 perl script, appears much like a native FORTRAN 77 compiler. If you are running Red Hat Linux, the following workaround is suggested (as superuser)

mv /usr/bin/f77 /usr/bin/f77.orig

ln -s /usr/bin/fort77 /usr/bin/f77

before compiling any acoustic propagation codes.

KRAKEN, by Michael Porter. Compiles with no modifications under Linux. [source] [manual]
OASES, by Henrik Schmidt. Also compiles without modification under Linux. [source] [manual]
MIT modified version of Michael Collins' RAM. Also compiles without modification under Linux. [Intel source] [Alpha source]
Ocean Acoustics Library at SAIC
How to get KRAKEN and RAM to agree in range-independent environments.
How to select the correct mesh parameters in RAM.

Propagation Code Wrappers

Most everyone has their own way of reading the output files which KRAKEN or other propagation codes generate. Many people use MATLAB functions to do this, but I personally find MATLAB to be cumbersome when processing large data sets. So, as part of my SM thesis (way back in 1995/1996), I wrote my own set of C function codes to read KRAKEN and PRUFER output files. The idea was to have a standard C interface to read the normal modes, and calculate Green's function, all using C.

Most of the code was developed using the Linux operating system, and is used primarily on Intel powered Linux boxes. The code was designed to read KRAKEN files generated by Compaq Alpha-based workstations runing Digital Unix, or generated on Intel Linux boxes.

Complex Math Library

The first, and most basic library is the complex math library. I know; who has not written their own complex math library in C? Still, I had to start somewhere and stupidly re-invented the wheel. This code forms the basis for all of the complex math in my code, so it is kind of important. It is composed of two files: the header file, and the program code.

Here is a description of each function. a and b are double precision floating point numbers, x, y, and z are complex numbers, and n is an integer:

y = make_complex(a, b) Takes real arguments a (real) and b (imaginary) and returns a complex number, y. (y = a + ib)

a = complex_real(y) Takes a complex argument, y and returns its real part, a. (a = Re(y))

b = complex_imag(y) Takes a complex argument, y, and returns its imaginary part, b. (b = Im(y))

z = complex_add(y, x) Adds complex argument y to complex argument x, and returns the result in z. (z = y + x)

z = complex_subtract(y, x) Subtracts complex argument y from complex argument x, and returns the result in z. (z = y - x)

m = complex_magnitude(y) Takes the magnitude of complex argument y and returns it in m. (m = abs(y))

p = complex_phase(y) Takes the phase (in radians) of complex argument y and returns it in p. (p = angle(y))

z = complex_mult (y, x) Multiplies complex arguments y and x together and returns the result in z. (z = y * x)

z = complex_div (y, x) Divides complex argument x into complex argument y, returning the result in z. (z = y / x)

z = complex_div_component (y, a) Divides real argument a into complex argument y, returning the result in z. (z = y / a)

z = complex_component_div (a, x) Divides complex argument x into real argument a, returning the result in z. (z = a / x)

z = complex_exp (y) Raises e to the complex argument y, and returns the result in z. (z = exp(y))

z = complex_sqrt (y) Takes the square root of complex argument y and returns the result in z. (z = sqrt(y))

z = complex_real_pow(y, a) Raises complex argument y to the real power a, and returns the result in z. (z = y^a)

z = complex_conj(y) Returns the complex conjugate of complex argument y in z. (z = y^*)

z = cvsum (*y, n) Sums the n-element complex array y and returns the result in z.

Complex Matrix Library

After the complex math library was written, the next wheel to re-invent was a complex matrix library. This provided for the standard matrix operators, as well as matrix inversion. The inversion routine was adapted from the real matrix inversion routine in Numerical Recipies. The matrix library requires the complex math library to work. The header file and source code are here (matrix_ops.c and nrzc.c), as well as a description for each function. To extract the eigenvalues and eigenvectors of a Hermitian matrix, one needs the files eig.c and eig_wrapper.c.

Before using a generic matrix function, one must allocate a matrix. For instance, if one wanted to add complex matrices A and B to result in matrix C, one would have to allocate matrices A, B, and C before executing the add_complex_matrix() function. This is true for all complex matrix functions listed here. Additionally, all matrix manipulation routines return a signed integer. If the return value is zero, this indicates successful execution of the requested operation. If the return value is -1, then some sort of error occurred, which is usually printed out on stderr.

a = allocate_complex_matrix (r, c) Allocates an r by c (row by column) complex matrix and returns it in A.

free_complex_matrix(a) Frees the memory allocated by complex matrix A, returning it to the operating system. No further references to A should be made unless it is re-allocated.

reshape_complex_matrix (a, r, c) Re-shapes complex matrix A into an r by c (row by column) complex matrix.

copy_complex_matrix (a, b) Copies the elements of complex matrix A into (previously allocated) complex matrix B.

transpose_complex_matrix (a) Transposes the elements of complex matrix A. (A = A^T)

conjugate_complex_matrix (a) Takes the complex conjugate of the elements of complex matrix A. (A = A^*)

hermitian_complex_matrix(a) Takes the complex conjugate of the elements of complex matrix A, and then transposes the complex matrix. (A = (A^*)^T)

add_complex_matrix (a, b, d) Adds complex matrices A and B, and stores the result in the elements of complex matrix D. (D = A + B)

subtract_complex_matrix (a, b, d) Subtracts complex matrix B from complex matrix A and stores the result in the elements of complex matrix D. (D = A - B)

scalar_multiply_complex_matrix (a, b, f) Multiplies the elements of complex matrix A by the scalar complex number f, and stores the result in complex matrix B. (B = f * A)

multiply_complex_matrix (a, b, d) Multiplies complex matrix A with complex matrix B and stores the result in complex matrix D. (D = A * B)

invert_complex_matrix (a, b) Takes the inverse of complex matrix A and stores the result in complex matrix B. (B = A^-1)

complex_eig (a, d, v) Obtains the eigenvalues and eigenvectors of Hermitian matrix A. Stores the real eigenvalues in a double precision array d. Stores the complex eigenvectors in an array of N x 1 complex matrices, V. You must allocate d using calloc() and and V using calloc() and allocate_complex_matrix() before calling this routine.

write_complex_matrix (fp, a) Writes complex matrix A into the file referenced by file pointer fp (returned by the fopen() C function)

a = read_complex_matrix (fp) Allocates and reads complex matrix A from the file referenced by file pointer fp (returned by the fopen() C function)

KRAKEN parsing functions

The next step was to write a program which would parse an output file generated by KRAKEN. This required creating a set of new structures, which are outlined in the header file mode.h. The actual routines for parsing the output file are outlined below. They rely on the complex math library presented above.
These routines will also need some glue code and another header file to read FORTRAN data,

m = read_all_kraken_modes (filename, rd, nrd, comp, swapbytes)

Reads all modes from a KRAKEN output file.

filename: The name of the KRAKEN output file, type char *
rd: Receiver depths where one would like to have sampled the modes, type float *
nrd: Number of receivers specified in the rd array
comp: The component of the mode which we are looking for. Typically, one wants the normal component, so this would be a character array with the first element as 'N'
swapbytes: 1 if bytes need to be swapped, 0 if not. If one runs KRAKEN on a machine with different byte ordering than the machine one runs this code, one will need a 1 here.
m: The mode information, (type Mode *), defined in the header file mode.h

mh = read_kraken_mode_header (filename, swpabytes)

Reads the header from a KRAKEN output file. This routine is called by read_all_kraken_modes; it provides all the information which read_all_kraken_modes does, except the actual mode shape. For example, this routine would be called if the user wanted to extract eigenvalue information from the KRAKEN output file.

filename: The name of the KRAKEN output file, type char *
swapbytes: 1 if bytes need to be swapped, 0 if not. If one runs KRAKEN on a machine with different byte ordering than the machine one runs this code, one will need a 1 here.
mh: The mode header information, (type Modeheader *), defined in the header file mode.h

Last updated: 000629
Comments/Questions: pmd@alum.mit.edu

`y = make_complex(a, b)`	Takes real arguments a (real) and b (imaginary) and returns a complex number, y. (y = a + ib)
`a = complex_real(y)`	Takes a complex argument, y and returns its real part, a. (a = Re(y))
`b = complex_imag(y)`	Takes a complex argument, y, and returns its imaginary part, b. (b = Im(y))
`z = complex_add(y, x)`	Adds complex argument y to complex argument x, and returns the result in z. (z = y + x)
`z = complex_subtract(y, x)`	Subtracts complex argument y from complex argument x, and returns the result in z. (z = y - x)
`m = complex_magnitude(y)`	Takes the magnitude of complex argument y and returns it in m. (m = abs(y))
`p = complex_phase(y)`	Takes the phase (in radians) of complex argument y and returns it in p. (p = angle(y))
`z = complex_mult (y, x)`	Multiplies complex arguments y and x together and returns the result in z. (z = y x)*
`z = complex_div (y, x)`	Divides complex argument x into complex argument y, returning the result in z. (z = y / x)
`z = complex_div_component (y, a)`	Divides real argument a into complex argument y, returning the result in z. (z = y / a)
`z = complex_component_div (a, x)`	Divides complex argument x into real argument a, returning the result in z. (z = a / x)
`z = complex_exp (y)`	Raises e to the complex argument y, and returns the result in z. (z = exp(y))
`z = complex_sqrt (y)`	Takes the square root of complex argument y and returns the result in z. (z = sqrt(y))
`z = complex_real_pow(y, a)`	Raises complex argument y to the real power a, and returns the result in z. (z = y^a)
`z = complex_conj(y)`	Returns the complex conjugate of complex argument y in z. (z = y^)*
`z = cvsum (*y, n)`	Sums the n-element complex array y and returns the result in z.

`a = allocate_complex_matrix (r, c)`	Allocates an r by c (row by column) complex matrix and returns it in A.
`free_complex_matrix(a)`	Frees the memory allocated by complex matrix A, returning it to the operating system. No further references to A should be made unless it is re-allocated.
`reshape_complex_matrix (a, r, c)`	Re-shapes complex matrix A into an r by c (row by column) complex matrix.
`copy_complex_matrix (a, b)`	Copies the elements of complex matrix A into (previously allocated) complex matrix B.
`transpose_complex_matrix (a)`	Transposes the elements of complex matrix A. (A = A^T)
`conjugate_complex_matrix (a)`	Takes the complex conjugate of the elements of complex matrix A. (A = A^)*
`hermitian_complex_matrix(a)`	Takes the complex conjugate of the elements of complex matrix A, and then transposes the complex matrix. (A = (A^)^T)*
`add_complex_matrix (a, b, d)`	Adds complex matrices A and B, and stores the result in the elements of complex matrix D. (D = A + B)
`subtract_complex_matrix (a, b, d)`	Subtracts complex matrix B from complex matrix A and stores the result in the elements of complex matrix D. (D = A - B)
`scalar_multiply_complex_matrix (a, b, f)`	Multiplies the elements of complex matrix A by the scalar complex number f, and stores the result in complex matrix B. (B = f A)*
`multiply_complex_matrix (a, b, d)`	Multiplies complex matrix A with complex matrix B and stores the result in complex matrix D. (D = A B)*
`invert_complex_matrix (a, b)`	Takes the inverse of complex matrix A and stores the result in complex matrix B. (B = A^-1)
`complex_eig (a, d, v)`	Obtains the eigenvalues and eigenvectors of Hermitian matrix A. Stores the real eigenvalues in a double precision array d. Stores the complex eigenvectors in an array of N x 1 complex matrices, V. You must allocate d using `calloc()` and and V using `calloc()` and `allocate_complex_matrix()` before calling this routine.
`write_complex_matrix (fp, a)`	Writes complex matrix A into the file referenced by file pointer `fp` (returned by the `fopen()` C function)
`a = read_complex_matrix (fp)`	Allocates and reads complex matrix A from the file referenced by file pointer `fp` (returned by the `fopen()` C function)