function [GROUPS] = swarming2d(M);

% W is the half-axis length for my field; L is the total length, for both
% x and y dimensions.  Spac is the spacing constant, allowing 'side' evenly sized 
% spaces along each axis.
W = sqrt(9*100*pi);
L = 2*W;

numquads = 36;
side = sqrt(numquads);

spac = (L/side);

% The cell which will hold the xy coords for each quad is initialized, Q_SET, along
% with a counter, 'rep'.
Q_SET = cell(numquads,1);
rep = 1;

% Now, an n x n  grid of squares (quads) will be defined, overlayed on the field.
for n = -W : spac : (W - spac)              % Loop to define lower x-coord.    
    Q_x = [n; (n + spac)];                  % This defines upper and lower x for 
    % the current rep.
    
    for m = -W : spac : (W - spac)          % Loop over y-coords.
        Q_y = [m; (m + spac)];              % Define upper/lower y's.
        
        Q = [Q_x(1,1) Q_y(1,1); ...
                Q_x(2,1) Q_y(1,1); ...
                Q_x(1,1) Q_y(2,1); ...
                Q_x(2,1) Q_y(2,1)];
        
        chop(Q,5);
        
        Q_SET{rep,1} = Q;           % Upper and lower x's and y's are stored in cell.
        
        rep = rep + 1;              % Move to next rep (index into cell).
        
    end                     % End m loop.
end                         % End n loop.


% The next routine finds the neighboring quads for each quad i.  This is a long and 
% ridiculous sequence, but it works (i.e. there is probably a better way).  Neighbors
% to quad i are stored in Q_NBS{i}.
Q_NBS = cell(numquads,1);           % Initialize Q_NBS.
for i = 1:numquads                  % Start i loop, over each cell.
    for j = 1:numquads              % J loop, over each possible neighbor.
        
        if i ~= j    
            SP = intersect(round(Q_SET{i}),round(Q_SET{j}),'rows');  
            
            if  size(SP,1) > 0 
                SC = [1];                     
            else
                SC = [];                      % In all other cases, SC will be empty.
            end
            
        else 
            SC = [];
        end
        
        if SC ~= []                     % If SC is not empty, add j to the neighbor 
            % list of i.           
            Q_NBS{i} = [Q_NBS{i} j];
        end        
    end                             % End j loop.
end                                 % End i loop.


% Now for the 'how many are there?' loop.
Q_NUMS = cell(numquads,1);

for q = 1:numquads
    
    xu = max(Q_SET{q}(:,1));
    xl = min(Q_SET{q}(:,1));
    yu = max(Q_SET{q}(:,2));
    yl = min(Q_SET{q}(:,2));
    
    XIS = find( M(:,2)> xl & M(:,2)< xu);
    YIS = find( M(:,3)> yl & M(:,3)< yu);
    
    XYS = intersect(XIS,YIS);
    
    Q_NUMS{q} = XYS;
    
end


% Now to check the counts and deal with grouping.
G_QUADS = [];
for cell = 1:numquads
    
    num = size(Q_NUMS{cell},1);
    
    if num == 3
        NBS = Q_NBS{cell};
        
        isol = 0;
        for ind = 1:size(NBS,2)
            nb = NBS(1,ind);
            nb_size = size(Q_NUMS{nb},1);
            
            if nb_size < 3
                nbs_count = 0;
            else
                nbs_count = 1;
            end
            
            isol = isol + nbs_count;
        end
        
        if isol ~= 0
            G_QUADS = [G_QUADS cell];
        end
        
        
    elseif num > 3
        G_QUADS = [G_QUADS cell];
    end
end
G_QUADS;   

% Now, the indicated quads need to be sorted into contiguous groups.
rep = 1;
while size(G_QUADS,2) ~= 0;
    
    qi = G_QUADS(1,1);
    
    QI_NBS = Q_NBS{qi};
    
    MATCH = intersect(QI_NBS,G_QUADS);
    
    if MATCH == []
        GROUPS(rep).members = qi;
        rep = rep + 1;
        killi = find(G_QUADS == qi);
        G_QUADS(killi) = [];
        
        
    else
        tmp_qi = qi;       
        TMP2 = [tmp_qi QI_NBS];
        TMP = intersect(TMP2,G_QUADS);
        GROUPS(rep).members = [];
        
        while size(TMP,2) ~= 0
            GROUPS(rep).members = [GROUPS(rep).members tmp_qi];
            iqi = find(TMP <= tmp_qi);          
            TMP(iqi) = [];
            
            if size(TMP,2) ~= 0
                tmp_qi = TMP(1,1);
                TMP2 = [TMP Q_NBS{tmp_qi}];
                TMP = intersect(TMP2,G_QUADS);
            end
        end
        
        GRP = GROUPS(rep).members;
        for kill = 1:size(GRP,2)
            get_ind = find(G_QUADS == GRP(1,kill));
            G_QUADS(get_ind) = [];
        end
        
        rep = rep+1;
        
    end 
end