wip adding bessel methods

Author: cassidymwagner

src/fluidsf/__init__.py

@@ -6,12 +6,14 @@
 from .bin_data import bin_data
 from .calculate_advection_2d import calculate_advection_2d
 from .calculate_advection_3d import calculate_advection_3d
+from .calculate_cascades_from_bessel_2d import calculate_cascades_from_bessel_2d
 from .calculate_separation_distances import calculate_separation_distances
 from .calculate_separation_distances_3d import calculate_separation_distances_3d
 from .calculate_sf_maps_2d import calculate_sf_maps_2d
 from .calculate_structure_function_1d import calculate_structure_function_1d
 from .calculate_structure_function_2d import calculate_structure_function_2d
 from .calculate_structure_function_3d import calculate_structure_function_3d
+from .calculate_wavenumbers_2d import calculate_wavenumbers_2d
 from .generate_sf_maps_2d import generate_sf_maps_2d
 from .generate_structure_functions_1d import generate_structure_functions_1d
 from .generate_structure_functions_2d import generate_structure_functions_2d
@@ -24,21 +26,23 @@
 __version__ = "0.2.0"
 
 __all__ = (
-    "generate_sf_maps_2d",
-    "generate_structure_functions_1d",
-    "generate_structure_functions_2d",
-    "generate_structure_functions_3d",
-    "calculate_sf_maps_2d",
-    "calculate_structure_function_1d",
-    "calculate_structure_function_2d",
-    "calculate_structure_function_3d",
+    "bin_data",
     "calculate_advection_2d",
     "calculate_advection_3d",
+    "calculate_cascades_from_bessel_2d",
     "calculate_separation_distances",
     "calculate_separation_distances_3d",
+    "calculate_sf_maps_2d",
+    "calculate_structure_function_1d",
+    "calculate_structure_function_2d",
+    "calculate_structure_function_3d",
+    "calculate_wavenumbers_2d",
+    "generate_sf_maps_2d",
+    "generate_structure_functions_1d",
+    "generate_structure_functions_2d",
+    "generate_structure_functions_3d",
     "shift_array_1d",
     "shift_array_2d",
     "shift_array_3d",
     "shift_array_xy",
-    "bin_data",
 )

src/fluidsf/calculate_cascades_from_bessel_2d.py

@@ -0,0 +1,45 @@
+from scipy.integrate import cumulative_trapezoid
+from scipy.special import jv
+import pandas as pd
+
+
+def calculate_cascades_from_bessel_2d(
+    kr, df, xdiffs, ydiffs, dx, dy, sf_keys=["SF_advection_velocity_x"]
+):
+    """
+    Calculate the cascades from the structure functions using Bessel functions.
+
+    Parameters
+    ----------
+    kr: numpy.ndarray
+        The kr wavevector.
+    df: pandas.DataFrame
+        The structure functions in a DataFrame.
+    xdiffs: numpy.ndarray
+        The x separation distances.
+    ydiffs: numpy.ndarray
+        The y separation distances.
+    dx: float
+        The x grid spacing.
+    dy: float
+        The y grid spacing.
+    sf_keys: list
+        The keys of the structure functions to use. Defaults to
+        ["SF_advection_velocity_x"].
+
+    Returns
+    -------
+    numpy.ndarray:
+        The cascades.
+    """
+    bessels1 = [jv(1, kr[i] * xdiffs) for i in range(kr.shape[0])]
+    bessels2 = [jv(2, kr[i] * xdiffs) for i in range(kr.shape[0])]
+    bessels3 = [jv(3, kr[i] * xdiffs) for i in range(kr.shape[0])]
+
+    for key in sf_keys:
+        df[key] = df[key].fillna(0)
+
+    bessel_energy_integral_asf_v = [
+        -0.5 * kr[kk] * cumulative_trapezoid(df[key] * bessels1[kk], dx=dx[kk])
+        for kk in range(1, len(kr) + 1)
+    ]

src/fluidsf/calculate_wavenumbers_2d.py

@@ -0,0 +1,41 @@
+import numpy as np
+
+
+def calculate_wavevector_2d(xdiffs, ydiffs):
+    """
+    Calculate the kr wavevector from the separation distances.
+    Assumes a uniform grid with no periodic boundary conditions.
+
+    Parameters
+    ----------
+    xdiffs: numpy.ndarray
+        The x separation distances.
+    ydiffs: numpy.ndarray
+        The y separation distances.
+
+    Returns
+    -------
+    numpy.ndarray:
+        The kr wavevector.
+    """
+    nx = len(xdiffs)
+    ny = len(ydiffs)
+    dx = abs(xdiffs[1] - xdiffs[0])
+    dy = abs(ydiffs[1] - ydiffs[0])
+
+    kx_int = 1 / dx
+    ky_int = 1 / dy
+
+    kx = 2 * np.py * np.arange(-nx / 2, nx / 2) * (kx_int / nx)
+    ky = 2 * np.py * np.arange(-ny / 2, ny / 2) * (ky_int / ny)
+
+    kx_max_mat = max(max(-kx[: nx // 2 + 1]), max(ky))
+
+    dkx = 2 * np.pi / max(xdiffs)
+    dky = 2 * np.pi / max(ydiffs)
+
+    dkr_mat = np.sqrt(dkx**2 + dky**2)
+
+    kr = np.arange(dkr_mat / 2, kx_max_mat + dkr_mat, dkr_mat)
+
+    return kr