function particles = particles(N)
%
%		Particles places a specified N particles into the 
%	plane at random positions and with random charges; it then 
%	calculates their position for ten timesteps.  An Euler method 
%	of integration is used.  Run it:
%
%               particles(3);
%
   rand('state',sum(100*clock));
   X = 2 * (rand(3, N) - 0.5);
   M = 2 * (rand(1, N) - 0.5)

   for i = 1:10,
      X = X + force(X, M) * 0.1
   end
return;

function force = force(X, M)
   F = zeros(size(X));
   for i = 1:length(M),
      for j = 1:length(M),
         if (i ~= j), 
            distance = norm(X(:,j) - X(:, i));
            force = M(j) / distance^2;
            for k = 1:size(F(:,i)),
               F(k, i) = F(k, i) + force * (X(k, j) - X(k, i))/distance;
            end
         end
      end
   end
   force = F;
return;