colide(v,cell_n,...

function [v,crmax,selxtra,col] = colide(v,cell_n,...
                 index,Xref,ncell,crmax,tau,selxtra,coeff)
% colide - Function to process collisions in cells
% Inputs
%    v       - Velocities of the particles
%    cell_n  - Number of particles in a cell
%    index   - Index for drawing from Xref list
%    Xref    - Cross-reference list for drawing particles
%    ncell   - Number of cells in the system
%    crmax   - Estimated maximum relative speed in a cell
%    tau     - Time step
%    selxtra - Extra selections carried over from last timestep
%    coeff   - Coefficient in computing number of selected pairs 
% Outputs
%    v       - Updated velocities of the particles
%    crmax   - Updated maximum relative speed
%    selxtra - Extra selections carried over to next timestep
%    col     - Total number of collisions processed
%
col = 0;          % Count number of collisions
%%%totsel = 0;
for jcell=1:ncell
 number = cell_n(jcell);
 if( number > 1 )  % If more than one particle is in the cell
  %% Determine number of candidate collision pairs 
  %% are to be selected in this cell
  select = coeff*number^2*crmax(jcell) + selxtra(jcell);
  nsel = floor(select);          % Number of pairs to be selected
  selxtra(jcell) = select-nsel;  % Carry over any left-over fraction
  crm = crmax(jcell);            % Current maximum relative speed
  for isel=1:nsel
    % Pick two particles at random out of this cell
    k = floor(rand(1)*number);
    kk = rem(ceil(k+rand(1)*(number-1)),number);
    ip1 = Xref(k+index(jcell));      % First particle
    ip2 = Xref(kk+index(jcell));     % Second particle
    cr = norm( v(ip1,:)-v(ip2,:) );  % Their relative speed 
    if( cr > crm )         % If relative speed larger than crm
      crm = cr;            % Reset crm to larger value
    end
    if( cr/crmax(jcell) > rand(1) )    % Accept or reject collision
      col = col+1;                     % Collision counter
      vcm = 0.5*(v(ip1,:) + v(ip2,:)); % Center of mass velocity
      cos_th = 1 - 2*rand(1);          % Cosine and sine of 
      sin_th = sqrt(1 - cos_th^2);     % collision angle theta
      phi = 2*pi*rand(1);              % Collision angle phi
      vrel(1) = cr*cos_th;             % Compute post-collision 
      vrel(2) = cr*sin_th*cos(phi);    % relative velocity
      vrel(3) = cr*sin_th*sin(phi);
      v(ip1,:) = vcm + 0.5*vrel;       % Update post-collision
      v(ip2,:) = vcm - 0.5*vrel;       % velocities
    end
  end
  crmax(jcell) = crm;     % Update max relative speed 
%%%  totsel = totsel + nsel;
 end
end	     
%%%fprintf('Fraction of selected used = %g\n',col/totsel);
return;