function [ f, pot, kin ] = compute_old ( np, nd, pos, vel, mass )
  f = zeros ( nd, np );
  pot = 0.0;
  for i = 1 : np
    for j = 1 : np
      if ( i ~= j )
        for k = 1 : nd
          rij(k) = pos(k,i) - pos(k,j);
        end
        d = 0.0;
        for k = 1 : nd
          d = d + rij(k)^2;
        end
        d = sqrt ( d );
        d2 = min ( d, pi / 2.0 );
        pot = pot + 0.5 * sin ( d2 ) * sin ( d2 );
        for k = 1 : nd
          f(k,i) = f(k,i) - rij(k) * sin ( 2.0 * d2 ) / d;
        end
      end
    end
  end
  kin = 0.0;
  for k = 1 : nd
    for j = 1 : np
      kin = kin + vel(k,j)^2;
    end
  end
  kin = 0.5 * mass * kin;
  return
end
