On Apr 14, 3:18 am, markus <markusg.p... at googlemail.com> wrote:
> Hi,
>
> I have a quite large list (length ~ 10,000) of coordinates (positions
> of particles from a simulation), e.g. data={{x1,y1,z1},
> {x2,y2,z2},...}, and I am looking for fast way to calculate the "pair
> correlation" or "radial distribution function", which tells me, for
> each value of r, how many particles are separated by a distance r (+-
> epsilon).
> Usually, this can be calculated by something like:
>
> Do[Do[histogr[[ Ceiling[Norm[ data[[i]]-data[[j]] ] / delta] ]], {j, i
> +1, Length[data]}], {i, 1, Length[data]-1}]
>
> where histogr is just a list whose index corresponds to the distance
> r.
>
> Unfortunately I have found that for lists of length >1,000, this way
> of calculating is very slow (it can take ~minutes), compared to the
> equivalent C code, which proceeds in a few seconds. I have also tried
> to "Compile" the function, but the speed does not increase...
>
> Does anybody know some fast way of calculating the pair distribution
> function in Mathematica?
This will do 12 000 particles in 55 sec on an old 1.6 GHz Mac G5.
The time scales roughly as n^2. Note that it wants the x, y, and z
vectors passed separately. It returns a list of frequency counts
for bins whose upper bounds are w, 2*w, ..., m@w, where m*w is at
least twice the distance from the centroid to the point farthest
from the centroid.
func = Compile[{{x,_Real,1},{y,_Real,1},{z,_Real,1},{w,_Real}},
Module[{n,m,f}, n = Length@x; m = Ceiling[2./w *
Sqrt@Max[(x-Tr@x/n)^2 + (y-Tr@y/n)^2 + (z-Tr@z/n)^2]];
f = Table[0,{m}]; Do[Scan[f[[#]]++&, Ceiling[Sqrt[
(Drop[x,i]-x[[i]])^2 + (Drop[y,i]-y[[i]])^2 +
(Drop[z,i]-z[[i]])^2]/w]], {i,n-1}]; f]]