Code covered by the BSD License

Highlights fromVector algebra for arrays of any size, with array expansion enabled

from Vector algebra for arrays of any size, with array expansion enabled by Paolo de Leva
Multiple dot, cross, and outer products, cross divisions, norms, normalizations, projections, etc.

ocd(c, b, idC, idB)
```function Aorth = ocd(c, b, idC, idB)
%OCD  Orthogonal cross division.
%   This function exploits the BAXFUN engine (MATLAB Central, file #23084),
%   which enables multiple cross divisions and array expansion (AX).
%
%   When B and C are 3-element vectors (e.g. 31, 13, or 113 arrays):
%
%       Aorth = OCD(C, B) returns the orthogonal cross division C /o B.
%
%   More generally, when B and C are arrays of any size containing one or
%   more 3-element vectors:
%
%       Aorth = OCD(C, B) is equivalent to Aorth = OCD(C, B, IDC, IDB),
%       where IDA and IDB are the first dimensions of A and B whose length
%       is 3.
%
%       Aorth = OCD(C, B, DIM) is equivalent to
%       Aorth = OCD(C, B, IDC, IDB), where IDC = IDB = DIM.
%
%       Aorth = OCD(C, B, IDC, IDB) returns the orthogonal cross divisions
%       between the vectors contained in C along dimension IDC and those
%       contained in B along dimension IDB. These vectors must have 3
%       elements. C and B are viewed as "block arrays". IDC and IDB are
%       referred to as their "internal dimensions" (IDs). For instance, a
%       362 array may be viewed as an array containing twelve 3-element
%       blocks. In this case, its size is denoted by (3)62, and its ID is
%       1. Since AX is enabled, C and B may have different size, and IDC
%       may not coincide with IDB (see MULTIPROD).
%
%       Input and output format:
%           Array     Block size     Internal dimension
%           ---------------------------------------------------------------
%           C         3  (1-D)       IDC
%           B         3  (1-D)       IDB
%           A         3  (1-D)       MAX(IDC, IDB)
%           ---------------------------------------------------------------
%           If SIZE(B)==SIZE(C) and IDB==IDC, then SIZE(A)=SIZE(B)=SIZE(C).
%
%   Two main kinds of cross divisions exist: indefinite (ICD) and definite
%   (DCD). See the help text of function CROSSDIV for details. The
%   orthogonal cross division (OCD) is the simplest form of DCD. It is
%   equal to the known and unique orthogonal term of an ICD:
%                      C /o B = Aorth = B x C / (B * B)                 (1)
%   When Aorth coincides with A (see example),
%                                C /o B = A                             (2)
%
%   Minimum-input OCDs
%       One of the components of B or C can be omitted, provided it is
%       relative to an axis non-orthogonal to C or B. Since B and C are by
%       definition orthogonal, the missing component can be computed by
%       solving the equation B * C = 0. To perform a minimum-input OCD, use
%       Aorth = CROSSDIV(TYPE, C, B, pi/2, IDC, IDB), where TYPE has one of
%       the following values (see CROSSDIV for details):
%           '+t-Bx'   (Unknown component of vector B)
%           '+t-By'
%           '+t-Bz'
%           '+t-Cx'   (Unknown component of vector C)
%           '+t-Cy'
%           '+t-Cz'
%
%   Warning:
%       By definition, in an OCD (C /o B) the divisor (B) must be the
%       second operand of a cross product (A x B = C). When you use as
%       divisor the first operand (A) of the cross product, rather than the
%       second, the OCD does not yield the expected result (Borth). Since
%       A x B = - B x A, it yields the opposite result:
%                              C /o A = - Borth                         (3)
%
%   Reference:
%       de Leva, P (2008). Anticrossproducts and cross divisions.
%       Journal of Biomechanics, 8, 1790-1800.
%
%   Examples:
%       V is the derivative of a vector R which rotates with an angular
%       velocity OMEGA. V has two components, parallel and orthogonal to R
%       (V = Vpar + Vorth). Since Vorth = OMEGA x R (Poisson's formula),
%       then
%                           OMEGA = OCD(Vorth, R)      [OMEGA = Vorth /o R]
%       Since R x Vorth = R x V, this is equivalent to
%                           OMEGA = OCD(V, R)              [OMEGA = V /o R]
%       valid in 3D even when the magnitude of R is not constant (Vpar~=0).
%
%       If C and B are both (3)62 arrays of vectors,
%       C = OCD(C, B) is  a (3)62 array  of vectors.
%
%       A single vector C is divided by thirty vectors contained in B:
%       If  C is ............. a (3)1   vector,
%       and B is ............. a 56(3) array of 30 vectors,
%       C = OCD(C, B, 1, 3) is a 56(3) array of 30 vectors.
%
%            PROJECTION, REJECTION, TESTOCD.

% \$ Version: 3.0 \$
% CODE      by:                 Paolo de Leva (IUSM, Rome, IT) 2009 Feb 2
% COMMENTS  by:                 Code author                    2009 Feb 26
% OUTPUT    tested by:          Code author                    2009 Feb 26
% -------------------------------------------------------------------------

% Allow 2 to 4 input arguments
error( nargchk(2, 4, nargin) );

% Setting IDC and/or IDB
switch nargin
case 2
idC = find(size(c)==3, 1, 'first'); % First dim. of length 3
idB = find(size(b)==3, 1, 'first');
if isempty(idC) || isempty(idB)
error('OCD:InvalidSize',...
'C and B must have at least one dimension of length 3.');
end
case 3
idB = idC;
end

% 1 - Cross product (B x C)
%     (This will issue an error if SIZE(B,IDB)~=3 or SIZE(C,IDC)~=3)
BxC = cross2(b, c, idB, idC);

% 2 - Dot product (B * B)
BB = dot(b, b, idB);

% 3 - Orthogonal term of cross division C // B.
%     (BAXFUN replicates B2 3 times along its singleton dimension IDB)
shiftB = max(0, idC-idB);
Aorth = baxfun(@rdivide, BxC, BB, 0, shiftB);

% NOTE: For vectors with null magnitude, the latter divison (by zero) will
% cause MATLAB to issue a warning. The respective normalized vector will be
% composed of NaNs.```