Calculating Particle Density with GDL/IDL in Gadget2

[Philerd]

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 I am am not very familiar with it and have ran into some problems.

My first problem is I am unsure as I am 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 don't know if I've written it totally wrong or if there is an other way to do it. The count function doesn't 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 doesn't 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 I am sure I am going about it wrong.

 

[lofgran]

Thank you for your question. It is great to see that you are using Gadget2 for your simulations and are interested in analyzing the density distribution of the particles. I am a scientist with expertise in computational astrophysics and I would be happy to assist you with your questions.

Firstly, I would like to address your concern about the efficiency of your method for calculating the density. You are correct in saying that checking every particle with every other particle is not an efficient approach, especially for large numbers of particles. However, there is a well-established method called the "particle-mesh" method which is commonly used in simulations to calculate the density distribution. This method involves dividing the simulation volume into a grid and assigning the particles to their nearest grid points. The density at each grid point is then calculated by summing the masses of the particles assigned to that grid point and dividing by the volume of the grid cell. This method is much more efficient as it avoids the need to check every particle with every other particle. I would recommend exploring this method for your simulations.

Secondly, regarding your code, I would like to suggest a few improvements. Firstly, you can use the "where" function to find the indices of the particles within a certain radius instead of using a loop. This would be more efficient and would also give you the indices of the particles, which you can then use to calculate the density. Secondly, you can use the "count" function to count the number of particles within a certain radius instead of using the "where" function. This would give you the desired result of a single number for the density at each particle position. I have included a modified version of your code below for your reference:

x = pos(0,*) ; Sets the x,y,z positions from the position array
y = pos(1,*)
z = pos(2,*)

r = 5000 ; Radius around particle to check for other particles.

indices = where(sqrt((x-pos(0,*))^2 + (y-pos(1,*))^2 + (z-pos(2,*))^2) < r, count) ; Finds indices of particles within radius r
density = count(indices) ; Calculates the number of particles within radius r at each particle position.

I hope this helps. Please feel free to ask any further questions or clarifications. Good luck with your simulations!