I had to find an efficient way to plot the density profiles of spherical systems of point particles of varying mass.  Here's a summary:

First note that the mass enclosed in a sphere of radius $r$ is

$$\displaystyle M(r)=4\pi\int_0^r\rho(r')r'^2 dr'$$

I'm trying to find $\rho(r)$, so I'd like to get rid of the integral sign. Thus I differentiate with respect to $r$

$$\displaystyle \frac{dM(r)}{dr}=4\pi \rho(r)r^2$$

then 

$$\displaystyle \rho(r)=\frac{1}{4\pi r^2}\frac{dM(r)}{dr}.$$

The task has been reduced to the problem of computing ${dM(r)}/{dr}$, which is rather straightforward.  If a fast sorting algorithm, such as quicksort, is used to sort the particles by distance from the center of the spherical distribution (i.e. radius), then a mass array M(i) can be quickly constructed such that each entry i is the mass contained between $r$ and $r+i\delta r$ for $0\leq r\leq R$, where $R=n\delta r$ and $n$ is the number of elements in the array. Then a finite differencing scheme can be used to compute the derivative ${dM(r)}/{dr}$ at each point. Finish by properly handling the case for which the derivative at i is vastly different from the derivatives at the surrounding points i-1 and i+1. This problem can be mitigated by using (at least) a central differencing scheme and choosing an appropriate value of $\delta r$; it should be no smaller than the resolution limit of the particle distribution in the region where the density is being computed/plotted.

There are, of course, several other ways of computing the density profile.  In particular, if the potential has already been computed everywhere, then the density distribution can be inferred from that. The method described above seemed to be the best method for my particular code since the quicksort algorithm typically scales as $\mathcal{O}(N\log N)$ and all of the other operations either scale as $\mathcal{O}(N)$ or $\mathcal{O}(n)$.

UPDATE:  Here's a plot of the density profile using only 10,000 particles and no smoothing (blue line) compared to the exact analytic density profile from which the distribution was built (red line).  Density is in arbitrary units and is scaled by the central density, so $\rho = \mbox{const} \rho(r)/\rho_c$.

When cloud-in-cell smoothing was added (green line below), the result was roughly the same as the un-smoothed version, but with less random variation and some systematic differences; the gradient in the particle number density caused a slight bias toward higher density at larger radii and lower central density. This is expected behavior because the particle distribution has been "blurred" in this case, which softens the density gradient.

This entry was posted on Thursday, February 24th, 2011 at 12:28 pm and is filed under regular update. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.