Repository: uncscode/particula

Issue:
Add g12 Radial Distribution Function
#### $g_{12}$: Radial Distribution Function

The radial distribution function is given by:

$$
g_{12} = \left(\frac{\eta^2 + r_c^2}{R^2 + r_c^2}\right)^{C_1/2}
$$

Where $C_1$ and $r_c$ are derived based on droplet and turbulence properties.

##### $C_1$: Calculation

$$
C_1 = \frac{y(St)}{\left(|\mathbf{g}| / (v_k / \tau_k)\right)^{f_3(R_\lambda)}}
$$

$$
y(St) = -0.1988 St^4 + 1.5275 St^3 - 4.2942 St^2 + 5.3406 St
$$

$$
f_3(R_\lambda) = 0.1886 \exp\left(\frac{20.306}{R_\lambda}\right)
$$

Where:

$$ St = max(St_1, St_2) $$

Since the fitting for $y(St)$ was done for a limited range of St in DNS,
it should be set to zero for large $St$ when the function $y(St)$ becomes negative.

##### $r_c$: Expression

$$
\left(\frac{r_c}{\eta}\right)^2 = |St_2 - St_1| F(a_o, R_\lambda)
$$

solving for $r_c$:

$$
r_c = \eta \sqrt{|St_2 - St_1| F(a_o, R_\lambda)}
$$

where:

$$
a_{Og} = a_o + \frac{\pi}{8} \left(\frac{|\mathbf{g}|}{v_k / \tau_k}\right)^2
$$

$$
F(a_{Og}, R_\lambda) = 20.115 \left(\frac{a_{Og}}{R_\lambda}\right)^{1/2}
$$

Discussion:
sourcery-ai[bot]: <!-- Generated by sourcery-ai[bot]: start review_guide -->

🧙 Sourcery is reviewing your pull request!

---

<details>
<summary>Tips and commands</summary>

#### Interacting with Sourcery

- **Trigger a new review:** Comment `@sourcery-ai review` on the pull request.
- **Continue discussions:** Reply directly to Sourcery's review comments.
- **Generate a GitHub issue from a review comment:** Ask Sourcery to create an
  issue from a review comment by replying to it.
- **Generate a pull request title:** Write `@sourcery-ai` anywhere in the pull
  request title to generate a title at any time.
- **Generate a pull request summary:** Write `@sourcery-ai summary` anywhere in
  the pull request body to generate a PR summary at any time. You can also use
  this command to specify where the summary should be inserted.

#### Customizing Your Experience

Access your [dashboard](https://app.sourcery.ai) to:
- Enable or disable review features such as the Sourcery-generated pull request
  summary, the reviewer's guide, and others.
- Change the review language.
- Add, remove or edit custom review instructions.
- Adjust other review settings.

#### Getting Help

- [Contact our support team](mailto:support@sourcery.ai) for questions or feedback.
- Visit our [documentation](https://docs.sourcery.ai) for detailed guides and information.
- Keep in touch with the Sourcery team by following us on [X/Twitter](https://x.com/SourceryAI), [LinkedIn](https://www.linkedin.com/company/sourcery-ai/) or [GitHub](https://github.com/sourcery-ai).

</details>

<!-- Generated by sourcery-ai[bot]: end review_guide -->
Gorkowski: @sourcery-ai review

Patch:
diff --git a/particula/dynamics/coagulation/turbulent_dns_kernel/g12_radial_distribution_ao2008.py b/particula/dynamics/coagulation/turbulent_dns_kernel/g12_radial_distribution_ao2008.py
new file mode 100644
index 000000000..b8921eb66
--- /dev/null
+++ b/particula/dynamics/coagulation/turbulent_dns_kernel/g12_radial_distribution_ao2008.py
@@ -0,0 +1,208 @@
+"""
+Calculate the radial distribution function g_{12} for particles in a
+turbulent flow.
+"""
+
+from numpy.typing import NDArray
+import numpy as np
+
+from particula.util.validate_inputs import validate_inputs
+from particula.util.constants import STANDARD_GRAVITY
+from particula.util.self_broadcast import (
+    get_pairwise_sum_matrix,
+    get_pairwise_max_matrix,
+    get_pairwise_diff_matrix,
+)
+
+
+@validate_inputs(
+    {
+        "particle_radius": "positive",
+        "stokes_number": "positive",
+        "kolmogorov_length_scale": "positive",
+        "reynolds_lambda": "positive",
+        "normalized_accel_variance": "positive",
+        "kolmogorov_velocity": "positive",
+        "kolmogorov_time": "positive",
+    }
+)
+def get_g12_radial_distribution_ao2008(
+    particle_radius: NDArray[np.float64],
+    stokes_number: NDArray[np.float64],
+    kolmogorov_length_scale: float,
+    reynolds_lambda: float,
+    normalized_accel_variance: float,
+    kolmogorov_velocity: float,
+    kolmogorov_time: float,
+) -> NDArray[np.float64]:
+    # pylint: disable=too-many-arguments, too-many-positional-arguments
+    """
+    Compute the radial distribution function g_{12}.
+
+    The radial distribution function describes the clustering of particles
+    in a turbulent flow and is given by:
+
+        g_{12} = ((η² + r_c²) / (R² + r_c²))^(C_1/2)
+
+    - η (kolmogorov_length_scale) : Kolmogorov length scale [m]
+    - R (collision_radius) : Collision radius (sum of particle radii) [m]
+    - Stokes_1, Stokes_2 : Stokes numbers of the two particles [-]
+    - R_λ (reynolds_lambda) : Taylor-microscale Reynolds number [-]
+    - a_o (normalized_accel_variance) : Normalized acceleration variance [-]
+    - v_k (kolmogorov_velocity) : Kolmogorov velocity scale [m/s]
+    - τ_k (kolmogorov_time) : Kolmogorov timescale [s]
+    - g (gravitational_acceleration) : Gravitational acceleration [m/s²]
+
+    Arguments:
+    ----------
+    - particle_radius : Array of particle radii [m]
+    - stokes_number : Array of Stokes numbers of particles [-]
+    - kolmogorov_length_scale : Kolmogorov length scale [m]
+    - reynolds_lambda : Taylor-microscale Reynolds number [-]
+    - normalized_accel_variance : Normalized acceleration variance [-]
+    - kolmogorov_velocity : Kolmogorov velocity scale [m/s]
+    - kolmogorov_time : Kolmogorov timescale [s]
+
+    Returns:
+    --------
+    - Radial distribution function g_{12} [-]
+
+    References:
+    -----------
+    - Ayala, O., Rosa, B., & Wang, L. P. (2008). Effects of turbulence on
+        the geometric collision rate of sedimenting droplets. Part 2.
+        Theory and parameterization. New Journal of Physics, 10.
+        https://doi.org/10.1088/1367-2630/10/7/075016
+    """
+    collision_radius = get_pairwise_sum_matrix(particle_radius)
+    stokes_diff_matrix = get_pairwise_diff_matrix(stokes_number)
+    stokes_max_matrix = get_pairwise_max_matrix(stokes_number)
+
+    c1 = _calculate_c1(
+        stokes_max_matrix,
+        reynolds_lambda,
+        kolmogorov_velocity,
+        kolmogorov_time,
+    )
+
+    rc = _calculate_rc(
+        stokes_diff_matrix,
+        kolmogorov_length_scale,
+        normalized_accel_variance,
+        reynolds_lambda,
+        kolmogorov_velocity,
+        kolmogorov_time,
+    )
+
+    return (
+        (kolmogorov_length_scale**2 + rc**2) / (collision_radius**2 + rc**2)
+    ) ** (c1 / 2)
+
+
+def _calculate_c1(
+    stokes_number: NDArray[np.float64],
+    reynolds_lambda: float,
+    kolmogorov_velocity: float,
+    kolmogorov_time: float,
+) -> NDArray[np.float64]:
+    """
+    Compute C_1 based on Stokes number and turbulence properties.
+
+    C_1 is given by:
+
+        C_1 = y(St) / (|g| / (v_k / τ_k))^f_3(R_λ)
+
+    - y(St) = -0.1988 St^4 + 1.5275 St^3 - 4.2942 St^2 + 5.3406 St
+    - f_3(R_λ) = 0.1886 * exp(20.306 / R_λ)
+    - |g| : Gravitational acceleration [m/s²]
+    - v_k : Kolmogorov velocity scale [m/s]
+    - τ_k : Kolmogorov timescale [s]
+    """
+    y_stokes = _compute_y_stokes(stokes_number)
+    f3_lambda = _compute_f3_lambda(reynolds_lambda)
+    gravity_term = np.abs(STANDARD_GRAVITY) / (
+        kolmogorov_velocity / kolmogorov_time
+    )
+
+    return y_stokes / (gravity_term**f3_lambda)
+
+
+def _compute_y_stokes(
+    stokes_number: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    """
+    Compute y(St), ensuring values remain non-negative.
+
+    y(St) = -0.1988 St^4 + 1.5275 St^3 - 4.2942 St^2 + 5.3406 St
+
+    Ensures y(St) ≥ 0 (if negative, sets to 0).
+    """
+    y_st = (
+        -0.1988 * stokes_number**4
+        + 1.5275 * stokes_number**3
+        - 4.2942 * stokes_number**2
+        + 5.3406 * stokes_number
+    )
+    return np.maximum(y_st, 0)  # Ensures y(St) is non-negative
+
+
+def _compute_f3_lambda(reynolds_lambda: float) -> float:
+    """Compute f_3(R_lambda), an empirical turbulence factor."""
+    return 0.1886 * np.exp(20.306 / reynolds_lambda)
+
+
+def _calculate_rc(
+    stokes_diff_matrix: NDArray[np.float64],
+    kolmogorov_length_scale: float,
+    normalized_accel_variance: float,
+    reynolds_lambda: float,
+    kolmogorov_velocity: float,
+    kolmogorov_time: float,
+) -> NDArray[np.float64]:
+    # pylint: disable=too-many-arguments, too-many-positional-arguments
+    """
+    Compute r_c, the turbulence-driven correction to the collision kernel.
+
+    The equation is:
+
+    (r_c / η)^2 = |St_2 - St_1| * F(a_Og, R_λ)
+    """
+    stokes_difference = np.abs(stokes_diff_matrix)
+    a_og = _compute_a_og(
+        normalized_accel_variance,
+        kolmogorov_velocity,
+        kolmogorov_time,
+    )
+    f_value = _compute_f(a_og, reynolds_lambda)
+
+    return kolmogorov_length_scale * np.sqrt(stokes_difference * f_value)
+
+
+def _compute_a_og(
+    normalized_accel_variance: float,
+    kolmogorov_velocity: float,
+    kolmogorov_time: float,
+) -> float:
+    """
+    Compute aOg, which accounts for the effect of gravity on
+    turbulence-driven clustering.
+
+    - a_Og = a_o + (π / 8) * (|g| / (v_k / τ_k))^2
+    """
+    gravity_term = np.abs(STANDARD_GRAVITY) / (
+        kolmogorov_velocity / kolmogorov_time
+    )
+    return normalized_accel_variance + (np.pi / 8) * gravity_term**2
+
+
+def _compute_f(
+    a_og: float,
+    reynolds_lambda: float,
+) -> float:
+    """
+    Compute F(aOg, R_lambda), an empirical scaling factor for
+    turbulence effects.
+
+    - F(a_Og, R_λ) = 20.115 * (a_Og / R_λ)^0.5
+    """
+    return 20.115 * (a_og / reynolds_lambda) ** 0.5
diff --git a/particula/util/self_broadcast.py b/particula/util/self_broadcast.py
new file mode 100644
index 000000000..9a78612cf
--- /dev/null
+++ b/particula/util/self_broadcast.py
@@ -0,0 +1,105 @@
+"""
+Self broadcasting operation to generate a pairwise sum matrix.
+"""
+
+import numpy as np
+from numpy.typing import NDArray
+
+
+def get_pairwise_sum_matrix(
+    input_array: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    """
+    Compute the pairwise sum matrix for a given 1D NumPy array.
+
+    This function performs a self-broadcasting operation to generate a matrix
+    where each element (i, j) is computed as:
+
+        output_matrix[i, j] = input_array[i] + input_array[j]
+
+    Arguments:
+    ----------
+        - input_array : A 1D NumPy array of values.
+
+    Returns:
+    --------
+        - A 2D NumPy array where each element is the sum of two elements
+          from the input array.
+
+    Example:
+    --------
+    ``` py title="Example"
+    input_array = np.array([1, 2, 3])
+    get_pairwise_sum_matrix(input_array)
+    array([[2, 3, 4],
+           [3, 4, 5],
+           [4, 5, 6]])
+    ```
+    """
+    return input_array[:, np.newaxis] + input_array[np.newaxis, :]
+
+
+def get_pairwise_diff_matrix(
+    input_array: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    """
+    Compute the pairwise difference matrix for a given 1D NumPy array.
+
+    This function performs a self-broadcasting operation to generate a matrix
+    where each element (i, j) is computed as:
+
+        output_matrix[i, j] = input_array[i] - input_array[j]
+
+    Arguments:
+    ----------
+        - input_array : A 1D NumPy array of values.
+
+    Returns:
+    --------
+        - A 2D NumPy array where each element is the difference of two
+          elements from the input array.
+
+    Example:
+    --------
+    ``` py title="Example"
+    input_array = np.array([1, 2, 3])
+    get_pairwise_diff_matrix(input_array)
+    array([[ 0, -1, -2],
+           [ 1,  0, -1],
+           [ 2,  1,  0]])
+    ```
+    """
+    return input_array[:, np.newaxis] - input_array[np.newaxis, :]
+
+
+def get_pairwise_max_matrix(
+    input_array: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    """
+    Compute the pairwise maximum matrix for a given 1D NumPy array.
+
+    This function performs a self-broadcasting operation to generate a matrix
+    where each element (i, j) is computed as:
+
+        output_matrix[i, j] = max(input_array[i], input_array[j])
+
+    Arguments:
+    ----------
+        - input_array : A 1D NumPy array of values.
+
+    Returns:
+    --------
+        - A 2D NumPy array where each element is the maximum of two
+          elements from the input array.
+
+    Example:
+    --------
+    ``` py title="Example"
+    input_array = np.array([1, 2, 3])
+    get_pairwise_max_matrix(input_array)
+    array([[1, 2, 3],
+           [2, 2, 3],
+           [3, 3, 3]])
+    ```
+    """
+    return np.maximum(input_array[:, np.newaxis], input_array[np.newaxis, :])


Test patch:
diff --git a/particula/dynamics/coagulation/turbulent_dns_kernel/tests/g12_radial_distribution_ao2008_test.py b/particula/dynamics/coagulation/turbulent_dns_kernel/tests/g12_radial_distribution_ao2008_test.py
new file mode 100644
index 000000000..873410945
--- /dev/null
+++ b/particula/dynamics/coagulation/turbulent_dns_kernel/tests/g12_radial_distribution_ao2008_test.py
@@ -0,0 +1,140 @@
+"""
+Test the radial distribution of the coagulation kernel g12 using the
+AO2008 DNS model.
+
+By calling the get_g12_radial_distribution_ao2008 function, we also
+test the private functions used in the calculation of the radial
+distribution function.
+"""
+
+import pytest
+import numpy as np
+
+from particula.dynamics.coagulation.turbulent_dns_kernel.g12_radial_distribution_ao2008 import (
+    get_g12_radial_distribution_ao2008,
+)
+
+
+def test_get_g12_radial_distribution_ao2008_scalar():
+    """
+    Test get_g12_radial_distribution_ao2008 with a small array input.
+    """
+    particle_radius = np.array(
+        [10e-6, 20e-6, 30e-6]
+    )  # Particle radii in meters
+    stokes_number = np.array([0.5, 0.7, 1.2])  # Stokes numbers
+    kolmogorov_length_scale = 1e-6  # m
+    reynolds_lambda = 100  # Taylor-microscale Reynolds number
+    normalized_accel_variance = 11  # Normalized acceleration variance
+    kolmogorov_velocity = 0.01  # m/s
+    kolmogorov_time = 0.005  # s
+
+    # Compute g12
+    g12_matrix = get_g12_radial_distribution_ao2008(
+        particle_radius,
+        stokes_number,
+        kolmogorov_length_scale,
+        reynolds_lambda,
+        normalized_accel_variance,
+        kolmogorov_velocity,
+        kolmogorov_time,
+    )
+
+    # Check the shape
+    assert g12_matrix.shape == (
+        3,
+        3,
+    ), f"Expected shape (3,3), but got {g12_matrix.shape}"
+
+    # Ensure values are non-negative
+    assert np.all(g12_matrix >= 0), "Expected all values to be non-negative"
+
+
+def test_get_g12_radial_distribution_ao2008_invalid_inputs():
+    """
+    Test that get_g12_radial_distribution_ao2008 raises validation errors
+    for invalid inputs.
+    """
+    particle_radius = np.array([10e-6, 20e-6, 30e-6])
+    stokes_number = np.array([0.5, 0.7, 1.2])
+    kolmogorov_length_scale = 1e-6
+    reynolds_lambda = 100
+    normalized_accel_variance = 11
+    kolmogorov_velocity = 0.01
+    kolmogorov_time = 0.005
+
+    with pytest.raises(ValueError):
+        get_g12_radial_distribution_ao2008(
+            -particle_radius,  # Negative values should raise an error
+            stokes_number,
+            kolmogorov_length_scale,
+            reynolds_lambda,
+            normalized_accel_variance,
+            kolmogorov_velocity,
+            kolmogorov_time,
+        )
+
+    with pytest.raises(ValueError):
+        get_g12_radial_distribution_ao2008(
+            particle_radius,
+            -stokes_number,  # should raise an error
+            kolmogorov_length_scale,
+            reynolds_lambda,
+            normalized_accel_variance,
+            kolmogorov_velocity,
+            kolmogorov_time,
+        )
+
+    with pytest.raises(ValueError):
+        get_g12_radial_distribution_ao2008(
+            particle_radius,
+            stokes_number,
+            -kolmogorov_length_scale,  # should raise an error
+            reynolds_lambda,
+            normalized_accel_variance,
+            kolmogorov_velocity,
+            kolmogorov_time,
+        )
+
+    with pytest.raises(ValueError):
+        get_g12_radial_distribution_ao2008(
+            particle_radius,
+            stokes_number,
+            kolmogorov_length_scale,
+            -reynolds_lambda,  # should raise an error
+            normalized_accel_variance,
+            kolmogorov_velocity,
+            kolmogorov_time,
+        )
+
+
+def test_get_g12_radial_distribution_ao2008_edge_cases():
+    """
+    Test get_g12_radial_distribution_ao2008 with edge cases, such as
+    very small or very large Stokes numbers.
+    """
+    particle_radius = np.array([10e-6, 20e-6, 30e-6])
+    stokes_number = np.array(
+        [0.000001, 0.01, 2000.0]
+    )  # Includes extreme Stokes numbers
+    kolmogorov_length_scale = 1e-6
+    reynolds_lambda = 100
+    normalized_accel_variance = 11
+    kolmogorov_velocity = 0.01
+    kolmogorov_time = 0.005
+
+    g12_matrix = get_g12_radial_distribution_ao2008(
+        particle_radius,
+        stokes_number,
+        kolmogorov_length_scale,
+        reynolds_lambda,
+        normalized_accel_variance,
+        kolmogorov_velocity,
+        kolmogorov_time,
+    )
+
+    assert g12_matrix.shape == (
+        3,
+        3,
+    ), f"Expected shape (3,3), but got {g12_matrix.shape}"
+    assert np.all(g12_matrix >= 0), "Expected all values to be non-negative"
diff --git a/particula/util/tests/self_broadcast_test.py b/particula/util/tests/self_broadcast_test.py
new file mode 100644
index 000000000..e40e4ea4c
--- /dev/null
+++ b/particula/util/tests/self_broadcast_test.py
@@ -0,0 +1,118 @@
+"""
+Test of the self_broadcast module in the util package.
+"""
+
+import numpy as np
+
+from particula.util.self_broadcast import (
+    get_pairwise_sum_matrix,
+    get_pairwise_diff_matrix,
+    get_pairwise_max_matrix,
+)
+
+
+def test_get_pairwise_sum_matrix():
+    """
+    Test get_pairwise_sum_matrix with a small array input.
+    """
+    input_array = np.array([1, 2, 3], dtype=np.float64)
+
+    expected_output = np.array(
+        [[2, 3, 4], [3, 4, 5], [4, 5, 6]], dtype=np.float64
+    )
+
+    result = get_pairwise_sum_matrix(input_array)
+
+    assert (
+        result.shape == expected_output.shape
+    ), f"Expected shape {expected_output.shape}, got {result.shape}"
+    assert np.allclose(
+        result, expected_output
+    ), f"Expected {expected_output}, but got {result}"
+
+
+def test_get_pairwise_diff_matrix():
+    """
+    Test get_pairwise_diff_matrix with a small array input.
+    """
+    input_array = np.array([1, 2, 3], dtype=np.float64)
+
+    expected_output = np.array(
+        [[0, -1, -2], [1, 0, -1], [2, 1, 0]], dtype=np.float64
+    )
+
+    result = get_pairwise_diff_matrix(input_array)
+
+    assert (
+        result.shape == expected_output.shape
+    ), f"Expected shape {expected_output.shape}, got {result.shape}"
+    assert np.allclose(
+        result, expected_output
+    ), f"Expected {expected_output}, but got {result}"
+
+
+def test_get_pairwise_max_matrix():
+    """
+    Test get_pairwise_max_matrix with a small array input.
+    """
+    input_array = np.array([1, 2, 3], dtype=np.float64)
+
+    expected_output = np.array(
+        [[1, 2, 3], [2, 2, 3], [3, 3, 3]], dtype=np.float64
+    )
+
+    result = get_pairwise_max_matrix(input_array)
+
+    assert (
+        result.shape == expected_output.shape
+    ), f"Expected shape {expected_output.shape}, got {result.shape}"
+    assert np.allclose(
+        result, expected_output
+    ), f"Expected {expected_output}, but got {result}"
+
+
+def test_pairwise_functions_edge_cases():
+    """
+    Test edge cases including empty arrays and single-element arrays.
+    """
+    empty_array = np.array([], dtype=np.float64)
+    single_element_array = np.array([5], dtype=np.float64)
+
+    # Test empty array
+    assert get_pairwise_sum_matrix(empty_array).size == 0
+    assert get_pairwise_diff_matrix(empty_array).size == 0
+    assert get_pairwise_max_matrix(empty_array).size == 0
+
+    # Test single-element array
+    assert np.allclose(
+        get_pairwise_sum_matrix(single_element_array), np.array([[10.0]])
+    )
+    assert np.allclose(
+        get_pairwise_diff_matrix(single_element_array), np.array([[0.0]])
+    )
+    assert np.allclose(
+        get_pairwise_max_matrix(single_element_array), np.array([[5.0]])
+    )
+
+
+def test_pairwise_functions_large_array():
+    """
+    Test that the functions work correctly for a larger array.
+    """
+    input_array = np.arange(10, dtype=np.float64)  # [0, 1, 2, ..., 9]
+
+    sum_matrix = get_pairwise_sum_matrix(input_array)
+    diff_matrix = get_pairwise_diff_matrix(input_array)
+    max_matrix = get_pairwise_max_matrix(input_array)
+
+    assert sum_matrix.shape == (10, 10)
+    assert diff_matrix.shape == (10, 10)
+    assert max_matrix.shape == (10, 10)
+
+    # Check specific expected values
+    assert sum_matrix[0, 9] == 9  # 0 + 9
+    assert diff_matrix[0, 9] == -9  # 0 - 9
+    assert max_matrix[0, 9] == 9  # max(0, 9)
+    assert sum_matrix[5, 5] == 10  # 5 + 5
+    assert diff_matrix[5, 5] == 0  # 5 - 5
+    assert max_matrix[5, 5] == 5  # max(5, 5)