What about changing this to something like:
rbin = floor(Int,r/rbins)
I cannot test it now, but I have the impression that you do not need to search there.
Do you know if it is the loop or the ball tree that is taking most of the time? I don’t know if NearestNeighbour
is the most appropriate package to compute what you want.