function [ f, pot, kin ] = compute ( np, nd, pos, vel, mass )

  f = zeros( nd, np );
  pot = 0.0;
  pi2 = pi / 2.0;
 
  for i = 1 : np
    Ri = pos - repmat ( pos( :, i ), 1, np );    % array of vectors to 'i'
    D = sqrt ( sum ( Ri.^2 ) );                  % array of distances
    Ri = Ri( :, ( D > 0.0 ) ); 
    D = D( D > 0.0 );                            % save only pos values
    D2 = D .* ( D <= pi2 ) + pi2 * ( D > pi2 );  % truncate the potential.
    pot = pot + 0.5 * sum ( sin ( D2 ).^2 );     % accumulate pot. energy
    f( :, i) = Ri * ( sin( 2*D2 ) ./ D );        % force on particle 'i'
  end
 
  kin = 0.5 * mass * sum ( diag ( vel' * vel) ); % kinetic energy
     
  return
end
