Nice workshop!

I see you're interested in this topic, so let me give you a IMO surprising information (if you don't know that already):
You've written a simulation where all N^2 pairwise interactions are calculated in, of course, O(N^2) time. However, there is a method that lets you calculate all interactions in linear time, i.e. with O(N) calculations, google for "Fast Multipole Method".