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

% function [xyze, xyz_mat] = beam_xyz(vel, 'angle', 30, 'orientation', 'down')
%
%      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 'orientation' argument pairs are optional;
%                defaults are as above.  orientation must be 'up'
%                or 'down'.
%      'direction', 'inverse' converts from xyze to beam
%                   'normal' (default) from beam to xyze

beam_angle = 30*pi/180;
orientation_up = 0; % default: down
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, 'orientation'),
      if strcmp(value, 'up'), orientation_up = 1;
      elseif ~strcmp(value, 'down'),
         eval('help beam_enu');
         error(['Invalid value of orientation: ', 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

cc = 1/(4*cos(beam_angle));
ss = 1/(2*sin(beam_angle));
xyz_mat = [  ss,  0, -ss,   0 ; ...
             0 , ss,   0, -ss ; ...
            -cc, -cc, -cc, -cc; ...
             cc, -cc,  cc, -cc ];

if orientation_up,
   xyz_mat([2 3],:) = -xyz_mat([2 3],:); % reverse y and z
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 sv(1) ~= 4,
      error('First dimension of vel must be 4 (for 4 beams).')
   end
   xyze = xyz_mat * vel;
end

