...

function [ave_n,ave_u,ave_T,nsamp]= ...
     sampler(x,v,npart,ncell,L,ave_n,ave_u,ave_T,nsamp)
% sampler - function to sample density and velocities
% Inputs
%    x     - Particle positions
%    v     - Particle velocities
%    npart - Number of particles
%    ncell - Number of cells
%    L     - System size
%    ave_n - Accumulated ave. number of particles in a cell
%    ave_u - Accumulated average fluid velocity in a cell
%    ave_T - Accumulated average temperature in a cell
%    nsamp - Number of samples taken
% Outputs
%    ave_n, ave_u, ave_T, nsamp - Updated values
%
jx=ceil(ncell*x/L);
sum_n = zeros(ncell,1);
sum_v = zeros(ncell,3);
sum_v2 = zeros(ncell,1);
for ipart=1:npart
  jcell = jx(ipart);  % Particle ipart is in cell jcell
  sum_n(jcell) = sum_n(jcell)+1;
  sum_v(jcell,:) = sum_v(jcell,:) + v(ipart,:);
  sum_v2(jcell) = sum_v2(jcell) + ...
              v(ipart,1)^2 + v(ipart,2)^2 + v(ipart,3)^2;
end
ave_n = ave_n + sum_n;
for i=1:3
  sum_v(:,i) = sum_v(:,i)./sum_n;
end
ave_u = ave_u + sum_v;
ave_T = ave_T + sum_v2./sum_n - ...
         (sum_v(:,1).^2 + sum_v(:,2).^2 + sum_v(:,3).^2);
nsamp = nsamp + 1;
return;