pyfracval.correlation
=====================

.. py:module:: pyfracval.correlation

.. autoapi-nested-parse::

   Density-density correlation function :math:`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
   (:cite:p:`Moran2019FracVAL` §4.2) therefore rests on :math:`f(r)`, whose
   log-log slope should approach :math:`D_f - d` with :math:`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 :math:`r`, a copy of the aggregate is
   displaced by :math:`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:

   .. math::
       f(r^{(k)}) = \frac{1}{n_{or}}\sum_{n_{or}} V_{int}^{(k)} / V_a

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

   Two spheres of radii :math:`r_1, r_2` whose centers are :math:`d` apart
   share the lens volume

   .. math::
       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 :math:`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 :math:`n_{or}=300` affordable at :math:`N=1024`.





Module Contents
---------------

.. py:function:: sphere_intersection_volume(r1, r2, d)

   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).

   :param r1: Broadcastable arrays of the two radii and the center separation.
   :type r1: np.ndarray
   :param r2: Broadcastable arrays of the two radii and the center separation.
   :type r2: np.ndarray
   :param d: Broadcastable arrays of the two radii and the center separation.
   :type d: np.ndarray

   :returns: Intersection volume for each triple.
   :rtype: np.ndarray


.. py:function:: density_correlation(coords, radii, n_orientations = 100, n_radii = 40, delta = 3.5, rng = None, densities = None)

   Estimate the density-density correlation function of an aggregate.

   :param coords: The aggregate geometry.
   :type coords: np.ndarray
   :param radii: The aggregate geometry.
   :type radii: np.ndarray
   :param n_orientations: Random displacement directions averaged per radius. The paper uses
                          300; 100 is usually enough to see the slope and is 3x cheaper.
   :type n_orientations: int
   :param n_radii: Number of geometrically-spaced radii.
   :type n_radii: int
   :param delta: Largest radius as a multiple of Rg.
   :type delta: float
   :param densities: Only used for the Rg that sets the radius range; f(r) itself is a
                     purely geometric (volume) quantity and is unaffected by density.
   :type densities: np.ndarray, optional

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


.. py:function:: fit_correlation_slope(result, fit_lo_over_rp = 2.0, fit_hi_over_rg = 1.0)

   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.
   :rtype: dict