Función de distribución radial

Ya no necesitas buscar más en otras páginas porque llegaste al lugar correcto, poseemos la respuesta que deseas sin complicarte.

Solución:

Mostraré un enfoque simple y rápido para calcular la función de correlación de pares (función de distribución radial) para un sistema 2D de partículas puntuales:

radialDistributionFunction2D[pts_?MatrixQ, boxLength_Real, nBins_: 350] :=
 Module[gr, r, binWidth = boxLength/(2 nBins), npts = [email protected], rho,
  rho = npts/boxLength^2; (* area number density *)
  r, gr = HistogramList[(*compute and bin the distances between points of interest*)
             Flatten @ DistanceMatrix @ pts, 0.005, boxLength/4., binWidth];
  r = MovingMedian[r, 2]; (* take center of each bin as r *)
  gr = gr/(2 Pi r rho binWidth npts); (* normaliza g(r) *)
  Transpose[r, gr] (* combine r and g(r) *)
 ]

Así es como se usa:

rdf = radialDistributionFunction2D[pts, 1023.];
ListLinePlot[rdf, PlotRange ->0, 150, All, Mesh -> 80]

Tenga en cuenta que obtiene la normalización correcta de forma gratuita. Esto tomó alrededor de 1.2 segundos en mi máquina. He restringido el rango de la trama para mostrar las características interesantes.

“Descargo de responsabilidad”: Este es únicamente un intento de limpiar / acelerar el código de OP tal como se presentó, ya que no he intentado encontrar un algoritmo alternativo para hacer el cálculo.

Primero algunos comentarios:

• Es posible que desee localizar todas las variables intermedias que utiliza en un Block.
• Cuando los índices más profundos de Part es All, no necesitamos especificarlos. Es decir, p[[All,2,All,All]] es simplemente equivalente a p[[All,2]].
• Table no es una construcción de bucle per se: Table produce un tabla de valores, así que a menos que hagamos algo con el resultado (p. ej. tab = Table[...]), deberíamos usar algo más. En el caso de OP, podrían ser reemplazados por Do, pero como veremos, podemos eliminarlos por completo.

El cuello de botella parece ser las llamadas a HistogramList, lo que sucede nCentral veces, 4397 en el ejemplo. Hay dos cosas que podemos hacer con esta parte del código: convertir el procedimiento g=0;Do[...;g=g+h] en un equivalente funcional usando Foldy reemplazando HistogramList con una variante casera y más rápida. Primero miramos el último punto.

La especificación del contenedor 0, maxShell - 1, 1 da contenedores unitarios $[0,1)$, $[1,2)$, etc., and since shell has integer values, HistogramList just tallies the values in shell that are less than maxShell - 1. Thus we can

1. Pick the shell values that are less than maxShell - 1
2. Tally a Sorted version of the result
3. Turn it into a SparseArray with default element 0 and correct length

Here is the code:

histList[shell_, maxShell_] : = SparseArray[
  Rule @@ Transpose[Tally[1 +
    Sort[Pick[shell, Negative[shell - maxShell + 1]]]maxShell - 1]

los g=0;Do[...;g=g+h] parte ahora se puede hacer funcional usando

gAcc[g_, i_] := Block[dist, shell,
  dist = p[[1, i + 1 ;; All]] - p[[1, i]], p[[2, i + 1 ;; All]] - p[[2, i]];
  shell = Floor[Sqrt[dist[[1]]^2 + dist[[2]]^2]/dr];
  g + histList[shell, maxShell]
]

y finalmente

g = Fold[gAcc, 0, Range[nCentral - 1]];

En mi máquina, esto es aproximadamente 4 veces más rápido que el código de OP.

Finalmente, la última parte que modifica g se puede vectorizar, por lo que no es necesario ningún bucle. Simplemente:

shells = Range[maxShell - 1];
areaShells = Pi*(((shells + 1.0)*dr)^2 - (shells*dr)^2);
g = g/(nCentral*areaShells*pointDensity);

image = [email protected]





"http://i.stack.imgur.com/czhuI.png"; 
pts = ComponentMeasurements[[email protected][image, BilateralFilter[image, 4, 1]], 
                           "Centroid"][[All, 2]];

lpts = [email protected];

(* the following is slow and not  really needed, we could take the center 
   to be at ImageDimensions/2 *)

nm = NMinimize[Tr[#.# & /@  Thread[Subtract[pts, x0 + Cos[ArcTan[-y0 + #[[2]],
            -x0 + #[[1]]]], y0 + Sin[ArcTan[-y0 + #[[2]], -x0 + #[[1]]]] & /@ 
                                     pts]]], x0, 400, 600, y0, 400, 600]

(*center*)
(*9.29992*10^8, x0 -> 504.879, y0 -> 505.181*)

ctr = x0, y0 /. nm[[2]]
centeredPts = Subtract[#, ctr] & /@ pts;
xyLimits = Transpose[-1/2 #, 1/2 # &@ctr];
ptsInZone = Select[pts, xyLimits[[1, 1]] < #[[1]] < xyLimits[[1, 2]] &&
                        xyLimits[[2, 1]] < #[[2]] < xyLimits[[2, 2]] &];

f = Nearest[pts];

g[pt_, rad_] := [email protected][pt, lpts, rad]
ran[n_] := Range[n, 250, n]
ptsAtAllRanges[pt_, n_] := g[pt, #] & /@ ran[n]

scale = 2;
pp = Differences[Length /@ ptsAtAllRanges[#, scale]] & /@ ptsInZone;
ListLinePlot[(Mean /@ [email protected])/[email protected][scale], PlotRange -> All]

La escala horizontal se divide por 2:

Sección de Reseñas y Valoraciones

Si estás de acuerdo, eres capaz de dejar un tutorial acerca de qué le añadirías a este enunciado.

---

---

Utiliza Nuestro Buscador