function theResult = bm2geo(theBeamData, theElevations, theAzimuths, ...
								theHeading, thePitch, theRoll, ...
								theOrientation, theBlankingDistance)
% bm2geo -- Convert ADCP beam-data to geographic components.
%  bm2geo(theBeamData, theElevations, theAzimuths, theHeading,
%   thePitch, theRoll, 'theOrientation', theBlankingDistance)
%   converts theBeamData to geographic coordinates, using the
%   beam configuration (theElevations, theAzimuths) and other
%   orientation information. Solutions with fewer than four
%   beams in the corresponding depth-cell are returned as NaNs.
%   All angles in degrees.  TheOrientation = {'down' | 'up'}.
%
% This code performs all the angular operations needed
%  to convert ADCP beam-coordinate data to geographic
%  coordinates.  It does not make adjustments for the
%  speed-of-sound or the ADCP frequency, both of which
%  can be done before or after the present operations.
%
% "bm2geo" calls the "bm2dir" and "bm2xyze" functions.
 
% Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 28-Apr-1998 11:49:00.

% Reference: ADCP Coordinate Transformation: Formulas and
%  Calculations (technical manual, 26 pages), RD Insruments,
%  1997.

if nargin < 1, help(mfilename), return, end
if nargin < 2, theElevations = [-70 -70 -70 -70]; end
if nargin < 3, theAzimuths = [270 90 0 180]; end
if nargin < 4, theHeading = 0; end
if nargin < 5, thePitch = 0; end
if nargin < 6, theRoll = 0; end
if nargin < 7, theOrientation = 'down'; end
if nargin < 8, theBlankingDistance = 0; end

% One elevation value given.

for j = length(theElevations)+1:4
	theElevations(j) = theElevations(j-1);
end

[m, n] = size(theBeamData);

% Depth-bins, scaled for 20 or 30 degree nominal beam-angle,
%  based on the z-components of the beam-direction matrix.
%  See reference page 8.

RCF = 180 ./ pi;

if mean(abs(theElevations)) > 65
	s = sin(70 ./ RCF);
else
	s = sin(60 ./ RCF);
end

theBeamDirections = bm2dir(theElevations, theAzimuths, ...
						theHeading, thePitch, theRoll, ...
						theOrientation);

depth_scale = (theBeamDirections(:, 3) ./ s).';
depth_bins = (1:m).' * abs(depth_scale);

% Interpolate the beam data in nearest-neighbor
%  fashion.  (The ADCP uses this scheme internally
%  when recording in earth-coordinates.)

result = zeros(size(theBeamData));
temp = zeros(1, n);
for j = 1:n
	theBeamData(:, j) = interp1(depth_bins(:, j), ...
							theBeamData(:, j), (1:m).', 'nearest');
end

% Transformation matrix.  See reference page 10.

theTransformation = bm2xyze(theElevations, theAzimuths, ...
						theHeading, thePitch, theRoll, ...
						theOrientation);

% Convert the beam-data (four-columns) to
%  earth-coordinates via post-multiplication
%  by the transpose of the transformation matrix.

result = theBeamData * theTransformation.';

if nargout > 0
	theResult = result;
else
	disp(result)
end
