sorter(x,npart,ncell,L)

function [cell_n, index, Xref] = sorter(x,npart,ncell,L)
% sorter - Function to sort particles into cells
% Inputs
%    x     - Positions of particles
%    npart - Number of particles in the system
%    ncell - Number of cells in the system
%    L     - System size
% Output
%    cell_n - Number of particles in a cell
%    index  - Index into Xref array
%    Xref   - Cross-reference array
[xsort, Xref] = sort(x);
jx = ceil(ncell*x/L);      % Particle i is in cell jx(i)
cell_n = zeros(ncell,1);
for ipart=1:npart
  jcell = jx(ipart);
  cell_n(jcell) = cell_n(jcell) + 1;
end
index = cumsum(cell_n);
index(2:ncell) = index(1:ncell-1) + 1;
index(1) = 1;
return;