function [w,d] = cancelMoment( u, d, cancelDrift )
%w = cancelMoment( u, d, cancelDrift )
%   d is a set of velocities of particles at positions u.
%   Calculate the angular velocity about the centroid of u that will cancel
%   out the moment of the velocities d.
%   If cancelDrift is true (default is false), then d will also be
%   translated to make its average zero.
%
%   The moment of velocities (assuming the centroid is zero) is sum_i u_i
%   cross d_i.

    if nargin < 3
        cancelDrift = false;
    end
    c = sum(u,1)/size(u,1);
    u = u - repmat(c,size(u,1),1);
    ucrossd = sum( cross(u,d,2), 1 );
    uu1 = sum( u(:,1).*u(:,1) );
    uu2 = sum( u(:,2).*u(:,2) );
    uu3 = sum( u(:,3).*u(:,3) );
    u1u2 = sum( u(:,1).*u(:,2) );
    u1u3 = sum( u(:,1).*u(:,3) );
    u2u3 = sum( u(:,2).*u(:,3) );
    u2u1 = sum( u(:,2).*u(:,1) );
    u3u1 = sum( u(:,3).*u(:,1) );
    u3u2 = sum( u(:,3).*u(:,2) );
    M = [ uu2+uu3, -u2u1, -u3u1; ...
          -u1u2, uu3+uu1, -u3u2; ...
          -u1u3, -u2u3, uu1+uu2 ];
    w = ucrossd / M;
    bigw = repmat( w, size(u,1), 1 );
    d = d - cross( bigw, u, 2 );
    if (length(cancelDrift)==1) && cancelDrift
        d = d - repmat( sum(d,1)/size(d,1), size(d,1), 1 );
    else
        for i=1:3
            if cancelDrift(i)
                d(:,i) = d(:,i) - sum(d(:,i))/size(d,1);
            end
        end
    end
    
    % test.
    % We should have u X (d - wXu) = 0
%     bigw = repmat( w, size(u,1), 1 );
%     test = sum( cross( u, d, 2 ), 1 )
%     uXwXu = sum( cross( cross( u, bigw, 2 ), u, 2 ), 1 )
%     wM = w*M
end
