pyfracval.correlation#

Density-density correlation function \(f(r)\) for fractal aggregates.

The radius of gyration is a single number, and an aggregate can match it while having quite the wrong internal structure. FracVAL’s own validation ([Morán et al., 2019] §4.2) therefore rests on \(f(r)\), whose log-log slope should approach \(D_f - d\) with \(d = 3\) over the scaling range. This module implements that estimator so generated aggregates can be checked against the metric the paper actually used - notably densified ones, which are reshaped after generation and so have the most to prove.

Method (paper Eqs. 14-15)#

For each of a set of distances \(r\), a copy of the aggregate is displaced by \(r\) in a random direction, the volume shared by the aggregate and its copy is computed analytically, and the result is averaged over orientations and normalized by the aggregate’s own volume:

\[f(r^{(k)}) = \frac{1}{n_{or}}\sum_{n_{or}} V_{int}^{(k)} / V_a\]

with radii spaced geometrically, \(r^{(k)} = R_{min}(R_{max}/R_{min})^{k/(n_{it}-1)}\), from \(R_{min} = r_{p,geo}/10\) to \(R_{max} = \delta R_g\) (\(\delta = 3.5\) is ample).

Two spheres of radii \(r_1, r_2\) whose centers are \(d\) apart share the lens volume

\[V = \frac{\pi (r_1 + r_2 - d)^2\,(d^2 + 2d r_2 - 3r_2^2 + 2d r_1 + 6 r_1 r_2 - 3 r_1^2)}{12 d}\]

which is exact, so no binning or kernel choice enters the estimate.

Cost#

Naively this is \(O(N^2)\) per orientation per radius. Since only pairs closer than the sum of two radii contribute anything, a k-d tree prunes it to the handful of genuinely overlapping pairs, which is what makes the paper’s \(n_{or}=300\) affordable at \(N=1024\).

Module Contents#

pyfracval.correlation.sphere_intersection_volume(r1, r2, d)[source]#

Exact shared volume of sphere pairs (vectorized).

Handles all three regimes: disjoint (0), one fully containing the other (the smaller sphere’s whole volume), and partial overlap (the lens formula).

Parameters:

  • r1 (np.ndarray) – Broadcastable arrays of the two radii and the center separation.

  • r2 (np.ndarray) – Broadcastable arrays of the two radii and the center separation.

  • d (np.ndarray) – Broadcastable arrays of the two radii and the center separation.

Returns:

Intersection volume for each triple.

Return type:

np.ndarray

pyfracval.correlation.density_correlation(coords, radii, n_orientations=100, n_radii=40, delta=3.5, rng=None, densities=None)[source]#

Estimate the density-density correlation function of an aggregate.

Parameters:

  • coords (np.ndarray) – The aggregate geometry.

  • radii (np.ndarray) – The aggregate geometry.

  • n_orientations (int) – Random displacement directions averaged per radius. The paper uses 300; 100 is usually enough to see the slope and is 3x cheaper.

  • n_radii (int) – Number of geometrically-spaced radii.

  • delta (float) – Largest radius as a multiple of Rg.

  • densities (np.ndarray, optional) – Only used for the Rg that sets the radius range; f(r) itself is a purely geometric (volume) quantity and is unaffected by density.

Returns:

r (radii), f (correlation values), r_over_rp (radii normalized by the geometric-mean primary radius), rg, and rp_geo.

Return type:

dict

pyfracval.correlation.fit_correlation_slope(result, fit_lo_over_rp=2.0, fit_hi_over_rg=1.0)[source]#

Fit the log-log slope of f(r), which should approach Df - 3.

The fit window matters and is not a free choice. Below roughly 2 r_p the curve reflects single-particle overlap rather than aggregate structure (the paper notes this explicitly), and beyond about Rg the finite size of the aggregate cuts the power law off. Only the range between them carries the fractal signal, and for small aggregates that range can be too short to fit meaningfully - which is itself worth reporting rather than hiding.

Returns:

slope, df_estimate (slope + 3), r_squared, n_points used, and the window actually used.

Return type:

dict