function c = cross(a,b,dim)
% Embedded MATLAB Library function.
%
% Limitations:
% If supplied, argument 3 (DIM) must be a constant.

% Copyright 2002-2004 The MathWorks, Inc.

eml_assert(nargin > 0, 'Input argument "A" is undefined.');
eml_assert(nargin > 1, 'Input argument "B" is undefined.');
% Assert:  same ndims AND (same dimensions OR (dim not supplied and a and 
% b are 3-element row or column vectors)).
eml_assert( length(size(a))==length(size(b)) && ( all(size(a)==size(b)) || ...
    ( nargin<3 && numel(a)==3 && numel(b)==3 && length(size(a))==2 ) ), ... 
    'A and B must be same size.');
eml_assert(nargin < 3 || eml_is_const(dim), ...
    'Input argument "DIM" must be a constant');
eml_assert(nargin < 3 || (dim >= 1 && dim <= length(size(a)) && size(a,dim) == 3), ...
    'A and B must be of length 3 in the dimension in which the cross product is taken.');
eml_assert(isfloat(a) && isfloat(b), ...
    'Unsupported input class.  Must be single or double.');

if isempty(a)
    c = a;
    return;
end

if numel(a) == 3 && length(size(a)) == 2 && ...
        (nargin < 3 || size(a,dim) == 3 && size(b,dim) == 3)
    %Quick processing for the usual case.
    c1 = a(2)*b(3) - a(3)*b(2);
    c2 = a(3)*b(1) - a(1)*b(3);
    c3 = a(1)*b(2) - a(2)*b(1);
    if (size(a,1) == 3 && size(b,1) == 3)
        c = [c1; c2; c3];
    else
        c = [c1, c2, c3];
    end
else
    aSize = size(a);
    if nargin < 3
        dim = findDim3(aSize);
        eml_assert(dim >= 1, ...
            'A and B must have at least one dimension of length 3.');
    end
    if (dim >= 2)
        stride = prod(aSize(1:dim-1));
    else
        stride = 1;
    end
    iNext = stride*3;
    if (dim >= size(aSize,2))
        nHigh = 1;
    else
        nHigh = 1 + iNext*(prod(aSize(dim+1:end)) - 1);
    end
    if isa(a,'single') || isa(b,'single')
        cls = 'single';
    else
        cls = 'double';
    end
    if isreal(a) && isreal(b)
        c = zeros(aSize,cls);
    else
        c = complex(zeros(aSize,cls));
    end
    for iStart = 1:iNext:nHigh
        for i1 = iStart:iStart+stride-1
            i2 = i1 + stride;
            i3 = i2 + stride;
            % Calculate cross product
            c(i1) = a(i2)*b(i3) - a(i3)*b(i2);
            c(i2) = a(i3)*b(i1) - a(i1)*b(i3);
            c(i3) = a(i1)*b(i2) - a(i2)*b(i1);
        end
    end
end


function k = findDim3(x)
%Find the index of the first value 3 in the vector x.
len = size(x,2);
k = 0;
for m = 1:len
    if x(m) == 3
        k = m;
        break;
    end
end