function [xyze, xyz_mat] = beam_xyz(vel,  varargin);

% function [xyze, xyz_mat] = beam_xyz(vel, 'angle', 30, ...
%                            'geometry', 'convex', 'direction', 'normal')
%
%      Transform velocities from RDI ADCP.
%
%      vel can be: beam velocity from adp.m,
%                  a vector of beam coordinates
%                  a matrix, the columns of which are beam coord vectors
%
%      'angle' and 'geometry' argument pairs are optional;
%                defaults are as above.
%
%      'direction', 'inverse' converts from xyze to beam
%                   'normal' (default) from beam to xyze
%
%  00/01/31 EF
%         modified matlab/ladp/beam_xyz.m, which was written for
%         Sontek data.
%         Based on "ADCP Coordinate Transformation", Copyright 1988 by
%         RD Instruments. pp 10-11
%
%       'orientation' is assumed to be 'down'.  'up' is not implimented
%

beam_angle = 30*pi/180;
convex = 1; % default: convex
xyze_to_beam = 0; % default: beam to xyze

for ii = 1:2:length(varargin)
   name = varargin{ii}; value = varargin{ii+1};
   if strcmp(name, 'angle'), beam_angle = value*pi/180;
   elseif strcmp(name, 'geometry'),
      if strcmp(value, 'concave'), convex = 0;
      elseif ~strcmp(value, 'convex'),
         eval('help beam_enu');
         error(['Invalid value of geometry: ', value]);
      end
   elseif strcmp(name, 'direction')
      if strcmp(value, 'inverse')
         xyze_to_beam = 1;
      elseif ~strcmp(value, 'normal')
         error('Invalid direction; must be normal or inverse');
      end
   else
      eval('help beam_enu')
      error(['Invalid argument name: ', name]);
   end
end

a = 1/(2*sin(beam_angle));
b = 1/(4*cos(beam_angle));
d = a/sqrt(2);

xyz_mat = [  a, -a,  0,  0; ...
             0,  0, -a,  a; ...
             b,  b,  b,  b; ...
             d,  d, -d, -d ];

if ~convex,
   xyz_mat([1, 2], :) = -xyz_mat([1, 2], :);
end

if xyze_to_beam,
   xyz_mat = inv(xyz_mat);
end

sv = size(vel);
ndims = length(sv);
if ndims > 3,
   error('Too many dimensions in the first input argument')
elseif ndims == 3
   xyze = zeros(sv);
   for ii = 1:sv(3)
      xyze(:,:,ii) = (xyz_mat * squeeze(vel(:,:,ii)).').';
   end
   return
else % has 2 dimensions even if it is empty
   if isempty(vel),
      xyze = [];
      return
   end
   if sv(1) == 1, vel = vel(:); end
   if size(vel, 1) ~= 4,
      error('First dimension of vel must be 4 (for 4 beams).')
   end
   xyze = xyz_mat * vel;
end

