I hope this is the right place to ask this. First off I am using Gadget2 to do simulations involving large numbers of particles and I would like to be able to calculate the density of each of the particles. The idea being when the density distribution is plotted it should be clear to see if any clumps have been formed or not. Anyway Gadget2 outputs the data into snapshot files at certain time intervals depending on what is set with the following information, particle positions, velocities, masses and any other information you might require. To read these files I need to use GDL / IDL, although I have used this language to write all the initial condition files Im am not very familiar with it and have ran into some problems. My first problem is im unsure as im going about the right method to calculate the density I am only interested in the position array of the particles and finding how many particles are in a certain radius. Which assuming the x and y position of a particle to be x and y and the x and y position of all other particles in the array as x1 and y1 I have written as: r1=sqrt((x-x1)^2)-((y-y1)^2) where r1 is the distance from the particle being looked at. But to calculate this I need to check every particle with every other particle in the array each time which is not very efficient. Is there another way to go about doing this, for example a transformation I could perform on the position array? Another problem is I have tried to write a program to do the above but its just not working right. So main part of the program in IDL / GDL x(*)=pos(0,0:N) ; Sets the x,y,z positions from the position array where y(*)=pos(1,0:N) ; Where N is the number of particles z(*)=pos(2,0:N) r = 5000 ; Radius around particle to check for other particles. for i=0,N do begin r1=sqrt((x(i)-pos(0,*))^2)-((y(i)-pos(1,*))^2) density=where(r1 < r ,count) endfor I dont know if Ive written it totally wrong or if there is an other way to do it. The count function doesnt do what I want it to do. I wanted it to simply count the number of particles it found resulting in a single number for the density array at position i to correspond to the particle at position i. However though it doesnt and saves the whole number from the position array of all the particles it finds to meet the condition, is there a proper function that will do what I want or a better way of writing it? Any help on this would be great as im sure im going about it wrong.