function theResult = bm2xyze(theElevations, theAzimuths, ...
								theHeading, thePitch, theRoll, ...
								theOrientation)
% bm2xyze -- ADCP beam to X-Y-Z-E transformation.
%  bm2xyze(theElevations, theAzimuths) returns the transformation
%   matrix that converts beam data (4 columns) to X-Y-Z-E data by
%   pre-multiplication, using RDI conventions for ADCP measurements.
%   The beam-directions point away from the transponders, whereas
%   positive velocities point toward them.  The error-vector
%   coefficients are scaled to be approximately 0.5.  All angles
%   in degrees.
%  bm2xyze(theElevations, theAzimuths, theHeading, thePitch, theRoll,
%   'theOrientation') uses the additional orientation information as
%   well.  All angles in degrees.  TheOrientation is {'down' | 'up'}.
%  bm2xyze (no arguments) demonstrates itself by returning the
%   transformation for the conventional arrangement of beams
%   pointed downwards 20 degrees from vertical, at azimuths
%   of [270 90 0 180], other angles = 0.


%%% START USGS BOILERPLATE -------------%
% Use of this program is described in:
%
% Acoustic Doppler Current Profiler Data Processing System Manual 
% Jessica M. Côté, Frances A. Hotchkiss, Marinna Martini, Charles R. Denham
% Revisions by: Andrée L. Ramsey, Stephen Ruane
% U.S. Geological Survey Open File Report 00-458 
% Check for later versions of this Open-File, it is a living document.
%
% Program written in Matlab v7.1.0 SP3
% Program updated in Matlab 7.2.0.232 (R2006a)
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%
%%% END USGS BOILERPLATE --------------

  
% 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 12-May-1998 14:47:33.

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

if nargin < 1
	help(mfilename)
	theElevations = [-70 -70 -70 -70]
	theAzimuths = [270 90 0 180]
	theHeading = 0
	thePitch = 0
	theRoll = 0
	theOrientation = 'down'
	theTransformation = bm2xyze(theElevations, theAzimuths, ...
							theHeading, thePitch, theRoll, ...
							theOrientation);
	if nargout > 0
		theResult = theTransformation;
	else
		theTransformation
	end
	return
end

% Default arguments.

if nargin < 2, theAzimuths = [270 90 0 180]; end
if nargin < 3, theHeading = 0; end
if nargin < 4, thePitch = 0; end
if nargin < 5, theRoll = 0; end
if nargin < 6, theOrientation = 0; end

% One elevation value given.

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

% From X-Y-Z to beam directions.

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

% From beam directions to X-Y-Z coordinates.
%  See reference page 10.

theInverse = theBeamDirections \ eye(4, 4);

% Error-vector.  See reference page 10.

theErrorVector = [theInverse(:, 1:3) \ theInverse(:, 4); -1];
theErrorVector = theErrorVector / norm(theErrorVector);
% just as a point of information, the error velocity result is inverted
% from the RDI winadcp output, this corrects the sign.  MM 3/30/10
% theErrorVector = -theErrorVector;

% Append.

theTransformation = [theInverse; theErrorVector.'];

% % MM 3/20/10 reviewing the rotations
% T = ones(4,4).*NaN;
% c = 1; % convex
% a = 1/(2*sin(20*pi/180)); % 20 deg. head example
% b = 1/(4*cos(20*pi/180)); % 20 deg. head example
% d = a/sqrt(2);
% T(1,:) = [c*a  -1*c*a   0      0];
% T(2,:) = [0     0      -1*c*a  c*a];
% T(3,:) = [b     b       b      b];
% T(4,:) = [d     d      -1*d    -1*d];
% 
% % lets see what the RDI formula for M looks like, pg. 13 of 2007 edition
% M = ones(3,3).*NaN;
% % need radians for MATLAB
% H = theHeading.*pi/180; % heading has been corrected for declination before this point
% P = thePitch.*pi/180; % pitch corrected for internal sensors above
% R = theRoll.*pi/180; % corrected already for orientation
% CH = cos(H); SH = sin(H);
% CP = cos(P); SP = sin(P);
% CR = cos(R); SR = sin(R);
% M(1,:) = [    CH.*CR+SH.*SP.*SR SH.*CP     CH.*SR+SH.*SP.*CR];
% M(2,:) = [-1.*SH.*CR+CH.*SP.*SR CH.*CP -1.*SH.*SR-CH.*SP.*CR];
% M(3,:) = [-1.*CP.*SR            SP         CP.*CR];
% 
% TM = M*T(1:3,1:3); % which will ignore beam 4

% Result for outward beam-directions, with positive
%  velocities directed inwards.

result = -theTransformation;

% Output.

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