From 4eecf621bfaefab3908e6d01b0a9fffc04c37501 Mon Sep 17 00:00:00 2001 From: Terrence Tricco Date: Wed, 27 Dec 2023 23:19:50 -0330 Subject: [PATCH 1/7] 1D surface density profile --- docs/api.rst | 11 ++++ sarracen/__init__.py | 1 + sarracen/disc/__init__.py | 1 + sarracen/disc/surface_density.py | 95 ++++++++++++++++++++++++++++++++ 4 files changed, 108 insertions(+) create mode 100644 sarracen/disc/__init__.py create mode 100644 sarracen/disc/surface_density.py diff --git a/docs/api.rst b/docs/api.rst index d19a61c..397c264 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -80,3 +80,14 @@ The default smoothing kernel is the cubic spline. Additional smoothing kernels a kernels.CubicSplineKernel kernels.QuarticSplineKernel kernels.QuinticSplineKernel + + +Disc +---- + +Accretion disc analysis routines are in the disc module. + +.. autosummary:: + :toctree: api/ + + disc.surface_density diff --git a/sarracen/__init__.py b/sarracen/__init__.py index 6689e77..2f50f04 100644 --- a/sarracen/__init__.py +++ b/sarracen/__init__.py @@ -8,5 +8,6 @@ from .interpolate import interpolate_2d, interpolate_2d_line, interpolate_3d_proj, interpolate_3d_cross from .render import render, streamlines, arrowplot +import sarracen.disc __version__ = "1.2.2" diff --git a/sarracen/disc/__init__.py b/sarracen/disc/__init__.py new file mode 100644 index 0000000..1a5644f --- /dev/null +++ b/sarracen/disc/__init__.py @@ -0,0 +1 @@ +from ..disc.surface_density import surface_density diff --git a/sarracen/disc/surface_density.py b/sarracen/disc/surface_density.py new file mode 100644 index 0000000..d6fe70d --- /dev/null +++ b/sarracen/disc/surface_density.py @@ -0,0 +1,95 @@ +import numpy as np +import pandas as pd +from ..sarracen_dataframe import SarracenDataFrame + + +def _get_mass(data: 'SarracenDataFrame'): + if data.mcol == None: + if 'mass' not in data.params: + raise KeyError("'mass' column does not exist in this SarracenDataFrame.") + return data.params['mass'] + + return data[data.mcol] + + +def surface_density(data: 'SarracenDataFrame', r_in: float = None, r_out: float = None, + bins: int = 300, geometry: str = 'cylindrical', origin: list = None, + retbins: bool = False) -> np.ndarray: + """ + Calculates the 1D azimuthally-averaged surface density profile. + + The surface density profile is computed by segmenting the particles into + radial bins (rings) and dividing the total mass contained within each bin + by the area of its respective ring. + + Parameters + ---------- + data : SarracenDataFrame + Particle data, in a SarracenDataFrame. + r_in : float, optional + Inner radius of the disc. Defaults to the minimum r value. + r_out : float, optional + Outer radius of the disc. Defaults to the maximum r value. + bins : int, optional + Defines the number of equal-width bins in the range [r_in, r_out]. + Default is 300. + geometry : str, optional + Coordinate system to use to calculate the particle radii. Can be + either *spherical* or *cylindrical*. Defaults to *cylindrical*. + origin : array-like, optional + The x, y and z position around which to compute radii. Defaults to + [0, 0, 0]. + retbins : bool, optional + Whether to return the bin edges or not. Defaults to False. + + Returns + ------- + array + A NumPy array of length bins containing the surface density profile. + array, optional + The location of the bin edges. Only returned if *retbins=True*. + + Raises + ------ + ValueError + If the *geometry* is not *cylindrical* or *spherical*. + + See Also + -------- + The surface density averaging procedure for SPH is described in section + 3.2.6 of Lodato & Price, MNRAS (2010), `doi:10.1111/j.1365-2966.2010.16526.x + `_. + """ + + if origin is None: + origin = [0, 0, 0] + + if geometry == 'spherical': + r = np.sqrt((data[data.xcol] - origin[0]) ** 2 + + (data[data.ycol] - origin[1]) ** 2 + + (data[data.zcol] - origin[2]) ** 2) + elif geometry == 'cylindrical': + r = np.sqrt((data[data.xcol] - origin[0]) ** 2 + + (data[data.ycol] - origin[1]) ** 2) + else: + raise ValueError("geometry should be either 'cylindrical' or 'spherical'") + + if r_in is None: + r_in = r.min() + if r_out is None: + r_out = r.max() + + bin_locations = np.linspace(r_in, r_out, bins+1) + areas = np.pi * (bin_locations[1:] ** 2 - bin_locations[:-1] ** 2) + rbins = pd.cut(r, bin_locations) + + mass = _get_mass(data) + if type(mass) is pd.Series: + sigma = mass.groupby(rbins).sum() + else: + sigma = data.groupby(rbins).count().iloc[:, 0] * mass + + if retbins: + return (sigma / areas).to_numpy(), bin_locations + else: + return (sigma / areas).to_numpy() \ No newline at end of file From 62b38f81022edf12e04f455c71899ba4eda028f4 Mon Sep 17 00:00:00 2001 From: Terrence Tricco Date: Thu, 28 Dec 2023 13:03:18 -0330 Subject: [PATCH 2/7] mass check uses isinstance instead of type --- sarracen/disc/surface_density.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sarracen/disc/surface_density.py b/sarracen/disc/surface_density.py index d6fe70d..59fba27 100644 --- a/sarracen/disc/surface_density.py +++ b/sarracen/disc/surface_density.py @@ -14,7 +14,7 @@ def _get_mass(data: 'SarracenDataFrame'): def surface_density(data: 'SarracenDataFrame', r_in: float = None, r_out: float = None, bins: int = 300, geometry: str = 'cylindrical', origin: list = None, - retbins: bool = False) -> np.ndarray: + retbins: bool = False): """ Calculates the 1D azimuthally-averaged surface density profile. @@ -84,7 +84,7 @@ def surface_density(data: 'SarracenDataFrame', r_in: float = None, r_out: float rbins = pd.cut(r, bin_locations) mass = _get_mass(data) - if type(mass) is pd.Series: + if isinstance(mass, pd.Series): sigma = mass.groupby(rbins).sum() else: sigma = data.groupby(rbins).count().iloc[:, 0] * mass From 4ec6637722e927860ebd2d184ce88b89feea78f2 Mon Sep 17 00:00:00 2001 From: Terrence Tricco Date: Thu, 28 Dec 2023 13:07:18 -0330 Subject: [PATCH 3/7] mass and origin surface density unit tests --- sarracen/tests/disc/__init__.py | 0 sarracen/tests/disc/test_surface_density.py | 58 +++++++++++++++++++++ 2 files changed, 58 insertions(+) create mode 100644 sarracen/tests/disc/__init__.py create mode 100644 sarracen/tests/disc/test_surface_density.py diff --git a/sarracen/tests/disc/__init__.py b/sarracen/tests/disc/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/sarracen/tests/disc/test_surface_density.py b/sarracen/tests/disc/test_surface_density.py new file mode 100644 index 0000000..3b483c7 --- /dev/null +++ b/sarracen/tests/disc/test_surface_density.py @@ -0,0 +1,58 @@ +import pandas as pd +import numpy as np +from numpy.testing import assert_array_equal, assert_allclose +from sarracen import SarracenDataFrame +from sarracen.disc import surface_density +import pytest + + +def test_mass_equivalency(): + """ Column mass result should equal global params mass result. """ + + # randomly place particles + rng = np.random.default_rng(seed=5) + x = rng.random(100) + y = rng.random(100) + z = rng.random(100) + mass = [3.2e-4] * 100 + + sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, 'mass': mass}) + sigma1 = surface_density(sdf) + + sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z}, + params={'mass': 3.2e-4}) + sigma2 = surface_density(sdf) + + assert_array_equal(sigma1, sigma2) + +@pytest.mark.parametrize("geometry", ['cylindrical', 'spherical']) +def test_origin(geometry): + """ Expect same profile regardless of origin. """ + + rng = np.random.default_rng(seed=5) + x = rng.random(100) + y = rng.random(100) + z = rng.random(100) + + sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z}, + params={'mass': 3.2e-4}) + sigma_zero = surface_density(sdf, origin=[0.0, 0.0, 0.0], + bins=30, geometry=geometry) + + sdf = SarracenDataFrame(data={'x': x - 2e-8, + 'y': y + 2e-8, + 'z': z + 3e-9}, + params={'mass': 3.2e-4}) + sigma_eps = surface_density(sdf, origin=[2e-8, -2e-8, 3e-9], + bins=30, geometry=geometry) + + sdf = SarracenDataFrame(data={'x': x + 5.6e4, + 'y': y - 8.7e3, + 'z': z + 5.4e6}, + params={'mass': 3.2e-4}) + sigma_large = surface_density(sdf, origin=[5.6e4, -8.7e3, 5.4e6], + bins=30, geometry=geometry) + + assert_allclose(sigma_zero, sigma_eps, rtol=2e-7) + assert_allclose(sigma_zero, sigma_large, rtol=2e-7) + From 6bcda7cf18a835541bd1ce87bcea93b495c27741 Mon Sep 17 00:00:00 2001 From: Terrence Tricco Date: Wed, 10 Jan 2024 00:22:47 -0330 Subject: [PATCH 4/7] disc scale height and angular momenta calculations --- sarracen/disc/__init__.py | 2 +- sarracen/disc/surface_density.py | 183 +++++++++++++++++--- sarracen/sarracen_dataframe.py | 65 ++++++- sarracen/tests/disc/test_angular_momenta.py | 62 +++++++ sarracen/tests/disc/test_surface_density.py | 28 ++- 5 files changed, 313 insertions(+), 27 deletions(-) create mode 100644 sarracen/tests/disc/test_angular_momenta.py diff --git a/sarracen/disc/__init__.py b/sarracen/disc/__init__.py index 1a5644f..a67f332 100644 --- a/sarracen/disc/__init__.py +++ b/sarracen/disc/__init__.py @@ -1 +1 @@ -from ..disc.surface_density import surface_density +from ..disc.surface_density import surface_density, angular_momenta, scale_height diff --git a/sarracen/disc/surface_density.py b/sarracen/disc/surface_density.py index 59fba27..5730f5e 100644 --- a/sarracen/disc/surface_density.py +++ b/sarracen/disc/surface_density.py @@ -1,5 +1,6 @@ import numpy as np import pandas as pd +import sys from ..sarracen_dataframe import SarracenDataFrame @@ -12,8 +13,72 @@ def _get_mass(data: 'SarracenDataFrame'): return data[data.mcol] -def surface_density(data: 'SarracenDataFrame', r_in: float = None, r_out: float = None, - bins: int = 300, geometry: str = 'cylindrical', origin: list = None, +def _bin_particles_by_radius(data: 'SarracenDataFrame', + r_in: float = None, + r_out: float = None, + bins: int = 300, + geometry: str = 'cylindrical', + origin: list = None): + """ + Utility function to bin particles in discrete intervals by radius. + + Parameters + ---------- + data: SarracenDataFrame + The particle dataset. + r_in : float, optional + Inner radius of the disc. Defaults to the minimum r value. + r_out : float, optional + Outer radius of the disc. Defaults to the maximum r value. + bins : int, optional + Defines the number of equal-width bins in the range [r_in, r_out]. + Default is 300. + geometry : str, optional + Coordinate system to use to calculate the particle radii. Can be + either *spherical* or *cylindrical*. Defaults to *cylindrical*. + origin : array-like, optional + The x, y and z position around which to compute radii. Defaults to + [0, 0, 0]. + + Returns + ------- + rbins: Series + The radial bin to which each particle belongs. + bin_locations: ndarray + Locations of the bin edges. + """ + + if origin is None: + origin = [0, 0, 0] + + if geometry == 'spherical': + r = np.sqrt((data[data.xcol] - origin[0]) ** 2 + + (data[data.ycol] - origin[1]) ** 2 + + (data[data.zcol] - origin[2]) ** 2) + elif geometry == 'cylindrical': + r = np.sqrt((data[data.xcol] - origin[0]) ** 2 + + (data[data.ycol] - origin[1]) ** 2) + else: + raise ValueError("geometry should be either 'cylindrical' or 'spherical'") + + # should we add epsilon here? + if r_in is None: + r_in = r.min() - sys.float_info.epsilon + if r_out is None: + r_out = r.max() + sys.float_info.epsilon + + bin_locations = np.linspace(r_in, r_out, bins+1) + rbins = pd.cut(r, bin_locations) + + return rbins, bin_locations + + +def surface_density(data: 'SarracenDataFrame', + r_in: float = None, + r_out: float = None, + bins: int = 300, + geometry: str = 'cylindrical', + origin: list = None, retbins: bool = False): """ Calculates the 1D azimuthally-averaged surface density profile. @@ -61,27 +126,10 @@ def surface_density(data: 'SarracenDataFrame', r_in: float = None, r_out: float `_. """ - if origin is None: - origin = [0, 0, 0] + rbins, bin_locations = _bin_particles_by_radius(data, r_in, r_out, bins, + geometry, origin) - if geometry == 'spherical': - r = np.sqrt((data[data.xcol] - origin[0]) ** 2 - + (data[data.ycol] - origin[1]) ** 2 - + (data[data.zcol] - origin[2]) ** 2) - elif geometry == 'cylindrical': - r = np.sqrt((data[data.xcol] - origin[0]) ** 2 - + (data[data.ycol] - origin[1]) ** 2) - else: - raise ValueError("geometry should be either 'cylindrical' or 'spherical'") - - if r_in is None: - r_in = r.min() - if r_out is None: - r_out = r.max() - - bin_locations = np.linspace(r_in, r_out, bins+1) areas = np.pi * (bin_locations[1:] ** 2 - bin_locations[:-1] ** 2) - rbins = pd.cut(r, bin_locations) mass = _get_mass(data) if isinstance(mass, pd.Series): @@ -92,4 +140,95 @@ def surface_density(data: 'SarracenDataFrame', r_in: float = None, r_out: float if retbins: return (sigma / areas).to_numpy(), bin_locations else: - return (sigma / areas).to_numpy() \ No newline at end of file + return (sigma / areas).to_numpy() + + +def _calc_angular_momenta(data: 'SarracenDataFrame', + rbins: pd.Series, + bin_locations: np.ndarray, + unit_vector: bool): + """ + Utility function to calculate angular momenta of the disc. + + Parameters + ---------- + data: SarracenDataFrame + Particle data, in a SarracenDataFrame. + rbins: Series + The radial bin to which each particle belongs. + bin_locations: ndarray + Locations of the bin edges. + unit_vector: bool + Whether to convert the angular momenta to unit vectors. + Default is True. + + Returns + ------- + Lx, Ly, Lz: Series + The x, y and z components of the angular momentum per bin. + """ + + mass = _get_mass(data) + + Lx = data[data.ycol] * data[data.vzcol] - data[data.zcol] * data[data.vycol] + Ly = data[data.zcol] * data[data.vxcol] - data[data.xcol] * data[data.vzcol] + Lz = data[data.xcol] * data[data.vycol] - data[data.ycol] * data[data.vxcol] + + if isinstance(mass, float): + Lx = (mass * Lx).groupby(rbins).sum() + Ly = (mass * Ly).groupby(rbins).sum() + Lz = (mass * Lz).groupby(rbins).sum() + else: + Lx = (data[data.mcol] * Lx).groupby(rbins).sum() + Ly = (data[data.mcol] * Ly).groupby(rbins).sum() + Lz = (data[data.mcol] * Lz).groupby(rbins).sum() + + if unit_vector: + Lmag = 1.0 / np.sqrt(Lx ** 2 + Ly ** 2 + Lz ** 2) + + Lx = Lx * Lmag + Ly = Ly * Lmag + Lz = Lz * Lmag + + return Lx, Ly, Lz + + +def angular_momenta(data: 'SarracenDataFrame', + r_in: float = None, + r_out: float = None, + bins: int = 300, + geometry: str = 'cylindrical', + origin: list = None, + retbins: bool = False, + unit_vector: bool = True): + + rbins, bin_locations = _bin_particles_by_radius(data, r_in, r_out, bins, + geometry, origin) + + Lx, Ly, Lz = _calc_angular_momenta(data, rbins, bin_locations, unit_vector) + Lx, Ly, Lz = Lx.to_numpy(), Ly.to_numpy(), Lz.to_numpy() + + if retbins: + return Lx, Ly, Lz, bin_locations + else: + return Lx, Ly, Lz + + +def scale_height(data: 'SarracenDataFrame', + r_in: float = None, + r_out: float = None, + bins: int = 300, + geometry: str = 'cylindrical', + origin: list = None, + retbins: bool = False): + rbins, bin_locations = _bin_particles_by_radius(data, r_in, r_out, bins, + geometry, origin) + + Lx, Ly, Lz = _calc_angular_momenta(data, rbins, bin_locations, + unit_vector=True) + + zdash = rbins.map(Lx).to_numpy() * data[data.xcol] \ + + rbins.map(Ly).to_numpy() * data[data.ycol] \ + + rbins.map(Lz).to_numpy() * data[data.zcol] + + return zdash.groupby(rbins).std().to_numpy() diff --git a/sarracen/sarracen_dataframe.py b/sarracen/sarracen_dataframe.py index 59ff3a8..1c36f04 100644 --- a/sarracen/sarracen_dataframe.py +++ b/sarracen/sarracen_dataframe.py @@ -36,7 +36,8 @@ class SarracenDataFrame(DataFrame): """ _internal_names = pd.DataFrame._internal_names + ['_xcol', '_ycol', '_zcol', - '_mcol', '_rhocol', '_hcol'] + '_mcol', '_rhocol', '_hcol', + '_vxcol', '_vycol', '_vzcol'] _internal_names_set = set(_internal_names) _metadata = ['_params', '_units', '_kernel'] @@ -106,7 +107,10 @@ def __init__(self, data=None, params=None, *args, **kwargs): self._units = None self.units = Series([np.nan for _ in range(len(self.columns))]) - self._xcol, self._ycol, self._zcol, self._hcol, self._mcol, self._rhocol = None, None, None, None, None, None + self._xcol, self._ycol, self._zcol = None, None, None + self._hcol, self._mcol, self._rhocol = None, None, None + self._vxcol, self._vycol, self._vzcol = None, None, None + self._identify_special_columns() self._kernel = CubicSplineKernel() @@ -167,6 +171,15 @@ def _identify_special_columns(self): elif 'density' in self.columns: self.rhocol = 'density' + # Look for the keyword 'rho' or 'density' in the data. + if 'vx' in self.columns: + self.vxcol = 'vx' + if 'vy' in self.columns: + self.vycol = 'vy' + if 'vz' in self.columns: + self.vzcol = 'vz' + + def create_mass_column(self): """ Create a new column 'm', copied from the 'massoftype' dataset parameter. @@ -473,6 +486,54 @@ def rhocol(self, new_col: str): if new_col in self or new_col is None: self._rhocol = new_col + + @property + def vxcol(self): + """ + str : Label of the column which contains the x-component of the velocity. + + If this is set to a column which does not exist in the dataset, the column + label will remain set to the old value. + """ + return self._vxcol + + @vxcol.setter + def vxcol(self, new_col: str): + if new_col in self or new_col is None: + self._vxcol = new_col + + @property + def vycol(self): + """ + str : Label of the column which contains the y-component of the velocity. + + If this is set to a column which does not exist in the dataset, the column + label will remain set to the old value. + """ + return self._vycol + + @vycol.setter + def vycol(self, new_col: str): + if new_col in self or new_col is None: + self._vycol = new_col + + + @property + def vzcol(self): + """ + str : Label of the column which contains the z-component of the velocity. + + If this is set to a column which does not exist in the dataset, the column + label will remain set to the old value. + """ + return self._vzcol + + @vzcol.setter + def vzcol(self, new_col: str): + if new_col in self or new_col is None: + self._vzcol = new_col + + @property def kernel(self): """ diff --git a/sarracen/tests/disc/test_angular_momenta.py b/sarracen/tests/disc/test_angular_momenta.py new file mode 100644 index 0000000..e7cc8fe --- /dev/null +++ b/sarracen/tests/disc/test_angular_momenta.py @@ -0,0 +1,62 @@ +import pandas as pd +import numpy as np +from numpy.testing import assert_array_equal, assert_allclose +from sarracen import SarracenDataFrame +from sarracen.disc import angular_momenta +import pytest + + +def test_mass_equivalency(): + """ Column mass result should equal global params mass result. """ + + # randomly place particles + rng = np.random.default_rng(seed=5) + x = rng.random(100) + y = rng.random(100) + z = rng.random(100) + vx = rng.random(100) + vy = rng.random(100) + vz = rng.random(100) + mass = [3.2e-4] * 100 + + sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, + 'vx': vx, 'vy': vy, 'vz': vz, + 'mass': mass}) + Lx1, Ly1, Lz1 = angular_momenta(sdf) + + sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, + 'vx': vx, 'vy': vy, 'vz': vz}, + params={'mass': 3.2e-4}) + Lx2, Ly2, Lz2 = angular_momenta(sdf) + + assert_array_equal(Lx1, Lx2) + assert_array_equal(Ly1, Ly2) + assert_array_equal(Lz1, Lz2) + + +def test_parts_vs_whole(): + """ Profiles should be the same for matching bins. """ + + # randomly place particles + rng = np.random.default_rng(seed=5) + x = rng.random(100) + y = rng.random(100) + z = rng.random(100) + vx = rng.random(100) + vy = rng.random(100) + vz = rng.random(100) + mass = [3.2e-4] * 100 + + sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, + 'vx': vx, 'vy': vy, 'vz': vz, + 'mass': mass}) + Lx_in, Ly_in, Lz_in = angular_momenta(sdf, r_in=0.0, r_out=0.5, bins=100) + Lx_out, Ly_out, Lz_out = angular_momenta(sdf, r_in=0.5, r_out=1.0, bins=100) + Lx_all, Ly_all, Lz_all = angular_momenta(sdf, r_in=0.0, r_out=1.0, bins=200) + + assert_array_equal(Lx_in, Lx_all[:100]) + assert_array_equal(Lx_out, Lx_all[100:]) + assert_array_equal(Ly_in, Ly_all[:100]) + assert_array_equal(Ly_out, Ly_all[100:]) + assert_array_equal(Lz_in, Lz_all[:100]) + assert_array_equal(Lz_out, Lz_all[100:]) \ No newline at end of file diff --git a/sarracen/tests/disc/test_surface_density.py b/sarracen/tests/disc/test_surface_density.py index 3b483c7..2ff02a8 100644 --- a/sarracen/tests/disc/test_surface_density.py +++ b/sarracen/tests/disc/test_surface_density.py @@ -25,6 +25,7 @@ def test_mass_equivalency(): assert_array_equal(sigma1, sigma2) + @pytest.mark.parametrize("geometry", ['cylindrical', 'spherical']) def test_origin(geometry): """ Expect same profile regardless of origin. """ @@ -53,6 +54,29 @@ def test_origin(geometry): sigma_large = surface_density(sdf, origin=[5.6e4, -8.7e3, 5.4e6], bins=30, geometry=geometry) - assert_allclose(sigma_zero, sigma_eps, rtol=2e-7) - assert_allclose(sigma_zero, sigma_large, rtol=2e-7) + assert_allclose(sigma_zero, sigma_eps, atol=2e-8, rtol=0.0) + assert_allclose(sigma_zero, sigma_large, atol=2e-8, rtol=0.0) + + +def test_parts_vs_whole(): + """ Profiles should be the same for matching bins. """ + + # randomly place particles + rng = np.random.default_rng(seed=5) + x = rng.random(100) + y = rng.random(100) + z = rng.random(100) + vx = rng.random(100) + vy = rng.random(100) + vz = rng.random(100) + mass = [3.2e-4] * 100 + + sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, + 'vx': vx, 'vy': vy, 'vz': vz, + 'mass': mass}) + sigma_in = surface_density(sdf, r_in=0.0, r_out=0.5, bins=100) + sigma_out = surface_density(sdf, r_in=0.5, r_out=1.0, bins=100) + sigma_all = surface_density(sdf, r_in=0.0, r_out=1.0, bins=200) + assert_allclose(sigma_in, sigma_all[:100], atol=1e-15, rtol=0.0) + assert_allclose(sigma_out, sigma_all[100:], atol=1e-15, rtol=0.0) From 087bee2d1639a353a501b6f4fff11123c17b0fe7 Mon Sep 17 00:00:00 2001 From: Terrence Tricco Date: Wed, 10 Jan 2024 09:35:00 -0330 Subject: [PATCH 5/7] calculate honH --- sarracen/disc/__init__.py | 2 +- sarracen/disc/surface_density.py | 239 +++++++++++++++--- ...ar_momenta.py => test_angular_momentum.py} | 12 +- 3 files changed, 212 insertions(+), 41 deletions(-) rename sarracen/tests/disc/{test_angular_momenta.py => test_angular_momentum.py} (82%) diff --git a/sarracen/disc/__init__.py b/sarracen/disc/__init__.py index a67f332..c35d96f 100644 --- a/sarracen/disc/__init__.py +++ b/sarracen/disc/__init__.py @@ -1 +1 @@ -from ..disc.surface_density import surface_density, angular_momenta, scale_height +from ..disc.surface_density import surface_density, angular_momentum, scale_height, honH diff --git a/sarracen/disc/surface_density.py b/sarracen/disc/surface_density.py index 5730f5e..cab6700 100644 --- a/sarracen/disc/surface_density.py +++ b/sarracen/disc/surface_density.py @@ -44,7 +44,7 @@ def _bin_particles_by_radius(data: 'SarracenDataFrame', ------- rbins: Series The radial bin to which each particle belongs. - bin_locations: ndarray + bin_edges: ndarray Locations of the bin edges. """ @@ -67,10 +67,10 @@ def _bin_particles_by_radius(data: 'SarracenDataFrame', if r_out is None: r_out = r.max() + sys.float_info.epsilon - bin_locations = np.linspace(r_in, r_out, bins+1) - rbins = pd.cut(r, bin_locations) + bin_edges = np.linspace(r_in, r_out, bins+1) + rbins = pd.cut(r, bin_edges) - return rbins, bin_locations + return rbins, bin_edges def surface_density(data: 'SarracenDataFrame', @@ -126,10 +126,10 @@ def surface_density(data: 'SarracenDataFrame', `_. """ - rbins, bin_locations = _bin_particles_by_radius(data, r_in, r_out, bins, + rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, geometry, origin) - areas = np.pi * (bin_locations[1:] ** 2 - bin_locations[:-1] ** 2) + areas = np.pi * (bin_edges[1:] ** 2 - bin_edges[:-1] ** 2) mass = _get_mass(data) if isinstance(mass, pd.Series): @@ -138,17 +138,16 @@ def surface_density(data: 'SarracenDataFrame', sigma = data.groupby(rbins).count().iloc[:, 0] * mass if retbins: - return (sigma / areas).to_numpy(), bin_locations + return (sigma / areas).to_numpy(), bin_edges else: return (sigma / areas).to_numpy() -def _calc_angular_momenta(data: 'SarracenDataFrame', - rbins: pd.Series, - bin_locations: np.ndarray, - unit_vector: bool): +def _calc_angular_momentum(data: 'SarracenDataFrame', + rbins: pd.Series, + unit_vector: bool): """ - Utility function to calculate angular momenta of the disc. + Utility function to calculate angular momentum of the disc. Parameters ---------- @@ -156,10 +155,8 @@ def _calc_angular_momenta(data: 'SarracenDataFrame', Particle data, in a SarracenDataFrame. rbins: Series The radial bin to which each particle belongs. - bin_locations: ndarray - Locations of the bin edges. unit_vector: bool - Whether to convert the angular momenta to unit vectors. + Whether to convert the angular momentum to unit vectors. Default is True. Returns @@ -193,27 +190,92 @@ def _calc_angular_momenta(data: 'SarracenDataFrame', return Lx, Ly, Lz -def angular_momenta(data: 'SarracenDataFrame', - r_in: float = None, - r_out: float = None, - bins: int = 300, - geometry: str = 'cylindrical', - origin: list = None, - retbins: bool = False, - unit_vector: bool = True): +def angular_momentum(data: 'SarracenDataFrame', + r_in: float = None, + r_out: float = None, + bins: int = 300, + geometry: str = 'cylindrical', + origin: list = None, + retbins: bool = False, + unit_vector: bool = True): + """ + Calculates the angular momentum profile of the disc. - rbins, bin_locations = _bin_particles_by_radius(data, r_in, r_out, bins, - geometry, origin) + The profile is computed by segmenting the particles into radial bins + (rings) and summing the angular momentum of the particles within each + bin. - Lx, Ly, Lz = _calc_angular_momenta(data, rbins, bin_locations, unit_vector) + Parameters + ---------- + data : SarracenDataFrame + Particle data, in a SarracenDataFrame. + r_in : float, optional + Inner radius of the disc. Defaults to the minimum r value. + r_out : float, optional + Outer radius of the disc. Defaults to the maximum r value. + bins : int, optional + Defines the number of equal-width bins in the range [r_in, r_out]. + Default is 300. + geometry : str, optional + Coordinate system to use to calculate the particle radii. Can be + either *spherical* or *cylindrical*. Defaults to *cylindrical*. + origin : array-like, optional + The x, y and z position around which to compute radii. Defaults to + [0, 0, 0]. + retbins : bool, optional + Whether to return the bin edges or not. Defaults to False. + + Returns + ------- + array + A NumPy array of length bins containing the angular momentum profile. + array, optional + The location of the bin edges. Only returned if *retbins=True*. + + Raises + ------ + ValueError + If the *geometry* is not *cylindrical* or *spherical*. + """ + + rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, + geometry, origin) + + Lx, Ly, Lz = _calc_angular_momentum(data, rbins, unit_vector) Lx, Ly, Lz = Lx.to_numpy(), Ly.to_numpy(), Lz.to_numpy() if retbins: - return Lx, Ly, Lz, bin_locations + return Lx, Ly, Lz, bin_edges else: return Lx, Ly, Lz +def _calc_scale_height(data: 'SarracenDataFrame', + rbins: pd.Series): + """ + Utility function to calculate the scale height of the disc. + + Parameters + ---------- + data: SarracenDataFrame + Particle data, in a SarracenDataFrame. + rbins: Series + The radial bin to which each particle belongs. + + Returns + ------- + H: Series + The scale height of the disc. + """ + Lx, Ly, Lz = _calc_angular_momentum(data, rbins, unit_vector=True) + + zdash = rbins.map(Lx).to_numpy() * data[data.xcol] \ + + rbins.map(Ly).to_numpy() * data[data.ycol] \ + + rbins.map(Lz).to_numpy() * data[data.zcol] + + return zdash.groupby(rbins).std() + + def scale_height(data: 'SarracenDataFrame', r_in: float = None, r_out: float = None, @@ -221,14 +283,123 @@ def scale_height(data: 'SarracenDataFrame', geometry: str = 'cylindrical', origin: list = None, retbins: bool = False): - rbins, bin_locations = _bin_particles_by_radius(data, r_in, r_out, bins, + """ + Calculates the scale height, H, of the disc. + + The scale height, H, is computed by segmenting the particles into radial + bins (rings) and calculating the angular momentum profile of the disc. + Each particle takes the dot product of its position vector with the + angular momentum vector of its corresponding bin. The standard deviation + of this result per bin yields the scale height profile of the disc. + + Parameters + ---------- + data : SarracenDataFrame + Particle data, in a SarracenDataFrame. + r_in : float, optional + Inner radius of the disc. Defaults to the minimum r value. + r_out : float, optional + Outer radius of the disc. Defaults to the maximum r value. + bins : int, optional + Defines the number of equal-width bins in the range [r_in, r_out]. + Default is 300. + geometry : str, optional + Coordinate system to use to calculate the particle radii. Can be + either *spherical* or *cylindrical*. Defaults to *cylindrical*. + origin : array-like, optional + The x, y and z position around which to compute radii. Defaults to + [0, 0, 0]. + retbins : bool, optional + Whether to return the bin edges or not. Defaults to False. + + Returns + ------- + array + A NumPy array of length bins scale height, H, profile. + array, optional + The location of the bin edges. Only returned if *retbins=True*. + + Raises + ------ + ValueError + If the *geometry* is not *cylindrical* or *spherical*. + + See Also + -------- + :func:`angular_momentum` : Calculate the angular momentum profile of a + disc. + """ + + rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, + geometry, origin) + + H = _calc_scale_height(data, rbins).to_numpy() + + if retbins: + return H, bin_edges + else: + return H + + +def honH(data: 'SarracenDataFrame', + r_in: float = None, + r_out: float = None, + bins: int = 300, + geometry: str = 'cylindrical', + origin: list = None, + retbins: bool = False): + """ + Calculates /H, the averaged smoothing length divided by the scale height. + + The profile is computed by segmenting the particles into radial bins + (rings). The average smoothing length in each bin is divided by the scale + height as calculated for that bin. + + Parameters + ---------- + data : SarracenDataFrame + Particle data, in a SarracenDataFrame. + r_in : float, optional + Inner radius of the disc. Defaults to the minimum r value. + r_out : float, optional + Outer radius of the disc. Defaults to the maximum r value. + bins : int, optional + Defines the number of equal-width bins in the range [r_in, r_out]. + Default is 300. + geometry : str, optional + Coordinate system to use to calculate the particle radii. Can be + either *spherical* or *cylindrical*. Defaults to *cylindrical*. + origin : array-like, optional + The x, y and z position around which to compute radii. Defaults to + [0, 0, 0]. + retbins : bool, optional + Whether to return the bin edges or not. Defaults to False. + + Returns + ------- + array + A NumPy array of length bins containing the /H profile. + array, optional + The location of the bin edges. Only returned if *retbins=True*. + + Raises + ------ + ValueError + If the *geometry* is not *cylindrical* or *spherical*. + + See Also + -------- + :func:`scale_height` : Calculate the scale height of a disc. + """ + + rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, geometry, origin) - Lx, Ly, Lz = _calc_angular_momenta(data, rbins, bin_locations, - unit_vector=True) + H = _calc_scale_height(data, rbins).to_numpy() - zdash = rbins.map(Lx).to_numpy() * data[data.xcol] \ - + rbins.map(Ly).to_numpy() * data[data.ycol] \ - + rbins.map(Lz).to_numpy() * data[data.zcol] + mean_h = data.groupby(rbins)[data.hcol].mean().to_numpy() - return zdash.groupby(rbins).std().to_numpy() + if retbins: + return mean_h / H, bin_edges + else: + return mean_h / H \ No newline at end of file diff --git a/sarracen/tests/disc/test_angular_momenta.py b/sarracen/tests/disc/test_angular_momentum.py similarity index 82% rename from sarracen/tests/disc/test_angular_momenta.py rename to sarracen/tests/disc/test_angular_momentum.py index e7cc8fe..6d3cd45 100644 --- a/sarracen/tests/disc/test_angular_momenta.py +++ b/sarracen/tests/disc/test_angular_momentum.py @@ -2,7 +2,7 @@ import numpy as np from numpy.testing import assert_array_equal, assert_allclose from sarracen import SarracenDataFrame -from sarracen.disc import angular_momenta +from sarracen.disc import angular_momentum import pytest @@ -22,12 +22,12 @@ def test_mass_equivalency(): sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, 'vx': vx, 'vy': vy, 'vz': vz, 'mass': mass}) - Lx1, Ly1, Lz1 = angular_momenta(sdf) + Lx1, Ly1, Lz1 = angular_momentum(sdf) sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, 'vx': vx, 'vy': vy, 'vz': vz}, params={'mass': 3.2e-4}) - Lx2, Ly2, Lz2 = angular_momenta(sdf) + Lx2, Ly2, Lz2 = angular_momentum(sdf) assert_array_equal(Lx1, Lx2) assert_array_equal(Ly1, Ly2) @@ -50,9 +50,9 @@ def test_parts_vs_whole(): sdf = SarracenDataFrame(data={'x': x, 'y': y, 'z': z, 'vx': vx, 'vy': vy, 'vz': vz, 'mass': mass}) - Lx_in, Ly_in, Lz_in = angular_momenta(sdf, r_in=0.0, r_out=0.5, bins=100) - Lx_out, Ly_out, Lz_out = angular_momenta(sdf, r_in=0.5, r_out=1.0, bins=100) - Lx_all, Ly_all, Lz_all = angular_momenta(sdf, r_in=0.0, r_out=1.0, bins=200) + Lx_in, Ly_in, Lz_in = angular_momentum(sdf, r_in=0.0, r_out=0.5, bins=100) + Lx_out, Ly_out, Lz_out = angular_momentum(sdf, r_in=0.5, r_out=1.0, bins=100) + Lx_all, Ly_all, Lz_all = angular_momentum(sdf, r_in=0.0, r_out=1.0, bins=200) assert_array_equal(Lx_in, Lx_all[:100]) assert_array_equal(Lx_out, Lx_all[100:]) From e4f93b6b61d586634e632bd952f2ddf28fd3d104 Mon Sep 17 00:00:00 2001 From: Terrence Tricco Date: Wed, 10 Jan 2024 09:57:07 -0330 Subject: [PATCH 6/7] angular momentum profile considers origin --- sarracen/disc/surface_density.py | 57 +++++++++++++++++++++----------- 1 file changed, 38 insertions(+), 19 deletions(-) diff --git a/sarracen/disc/surface_density.py b/sarracen/disc/surface_density.py index cab6700..6ccbbfb 100644 --- a/sarracen/disc/surface_density.py +++ b/sarracen/disc/surface_density.py @@ -13,6 +13,13 @@ def _get_mass(data: 'SarracenDataFrame'): return data[data.mcol] +def _get_origin(origin: list) -> list: + if origin is None: + return [0.0, 0.0, 0.0] + else: + return origin + + def _bin_particles_by_radius(data: 'SarracenDataFrame', r_in: float = None, r_out: float = None, @@ -48,9 +55,6 @@ def _bin_particles_by_radius(data: 'SarracenDataFrame', Locations of the bin edges. """ - if origin is None: - origin = [0, 0, 0] - if geometry == 'spherical': r = np.sqrt((data[data.xcol] - origin[0]) ** 2 + (data[data.ycol] - origin[1]) ** 2 @@ -126,8 +130,9 @@ def surface_density(data: 'SarracenDataFrame', `_. """ + origin = _get_origin(origin) rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, - geometry, origin) + geometry, origin) areas = np.pi * (bin_edges[1:] ** 2 - bin_edges[:-1] ** 2) @@ -145,6 +150,7 @@ def surface_density(data: 'SarracenDataFrame', def _calc_angular_momentum(data: 'SarracenDataFrame', rbins: pd.Series, + origin: list, unit_vector: bool): """ Utility function to calculate angular momentum of the disc. @@ -155,6 +161,8 @@ def _calc_angular_momentum(data: 'SarracenDataFrame', Particle data, in a SarracenDataFrame. rbins: Series The radial bin to which each particle belongs. + origin: list + unit_vector: bool Whether to convert the angular momentum to unit vectors. Default is True. @@ -167,9 +175,13 @@ def _calc_angular_momentum(data: 'SarracenDataFrame', mass = _get_mass(data) - Lx = data[data.ycol] * data[data.vzcol] - data[data.zcol] * data[data.vycol] - Ly = data[data.zcol] * data[data.vxcol] - data[data.xcol] * data[data.vzcol] - Lz = data[data.xcol] * data[data.vycol] - data[data.ycol] * data[data.vxcol] + x_data = data[data.xcol].to_numpy() - origin[0] + y_data = data[data.ycol].to_numpy() - origin[1] + z_data = data[data.zcol].to_numpy() - origin[2] + + Lx = y_data * data[data.vzcol] - z_data * data[data.vycol] + Ly = z_data * data[data.vxcol] - x_data * data[data.vzcol] + Lz = x_data * data[data.vycol] - y_data * data[data.vxcol] if isinstance(mass, float): Lx = (mass * Lx).groupby(rbins).sum() @@ -238,10 +250,11 @@ def angular_momentum(data: 'SarracenDataFrame', If the *geometry* is not *cylindrical* or *spherical*. """ + origin = _get_origin(origin) rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, geometry, origin) - Lx, Ly, Lz = _calc_angular_momentum(data, rbins, unit_vector) + Lx, Ly, Lz = _calc_angular_momentum(data, rbins, origin, unit_vector) Lx, Ly, Lz = Lx.to_numpy(), Ly.to_numpy(), Lz.to_numpy() if retbins: @@ -251,7 +264,8 @@ def angular_momentum(data: 'SarracenDataFrame', def _calc_scale_height(data: 'SarracenDataFrame', - rbins: pd.Series): + rbins: pd.Series, + origin: list = None): """ Utility function to calculate the scale height of the disc. @@ -261,13 +275,16 @@ def _calc_scale_height(data: 'SarracenDataFrame', Particle data, in a SarracenDataFrame. rbins: Series The radial bin to which each particle belongs. + origin : array-like, optional + The x, y and z position around which to compute radii. Defaults to + [0, 0, 0]. Returns ------- H: Series The scale height of the disc. """ - Lx, Ly, Lz = _calc_angular_momentum(data, rbins, unit_vector=True) + Lx, Ly, Lz = _calc_angular_momentum(data, rbins, origin, unit_vector=True) zdash = rbins.map(Lx).to_numpy() * data[data.xcol] \ + rbins.map(Ly).to_numpy() * data[data.ycol] \ @@ -330,10 +347,11 @@ def scale_height(data: 'SarracenDataFrame', disc. """ + origin = _get_origin(origin) rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, geometry, origin) - H = _calc_scale_height(data, rbins).to_numpy() + H = _calc_scale_height(data, rbins, origin).to_numpy() if retbins: return H, bin_edges @@ -342,12 +360,12 @@ def scale_height(data: 'SarracenDataFrame', def honH(data: 'SarracenDataFrame', - r_in: float = None, - r_out: float = None, - bins: int = 300, - geometry: str = 'cylindrical', - origin: list = None, - retbins: bool = False): + r_in: float = None, + r_out: float = None, + bins: int = 300, + geometry: str = 'cylindrical', + origin: list = None, + retbins: bool = False): """ Calculates /H, the averaged smoothing length divided by the scale height. @@ -392,10 +410,11 @@ def honH(data: 'SarracenDataFrame', :func:`scale_height` : Calculate the scale height of a disc. """ + origin = _get_origin(origin) rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, - geometry, origin) + geometry, origin) - H = _calc_scale_height(data, rbins).to_numpy() + H = _calc_scale_height(data, rbins, origin).to_numpy() mean_h = data.groupby(rbins)[data.hcol].mean().to_numpy() From d353ba561bcd0c9d2e60bd6a3365653bd9ecb445 Mon Sep 17 00:00:00 2001 From: Terrence Tricco Date: Wed, 10 Jan 2024 12:59:20 -0330 Subject: [PATCH 7/7] retbins returns the midpoints of each bin, instead of bin edges --- sarracen/disc/surface_density.py | 56 +++++++++++++++++++------------- 1 file changed, 33 insertions(+), 23 deletions(-) diff --git a/sarracen/disc/surface_density.py b/sarracen/disc/surface_density.py index 6ccbbfb..735146e 100644 --- a/sarracen/disc/surface_density.py +++ b/sarracen/disc/surface_density.py @@ -44,7 +44,7 @@ def _bin_particles_by_radius(data: 'SarracenDataFrame', Coordinate system to use to calculate the particle radii. Can be either *spherical* or *cylindrical*. Defaults to *cylindrical*. origin : array-like, optional - The x, y and z position around which to compute radii. Defaults to + The x, y and z centre point around which to compute radii. Defaults to [0, 0, 0]. Returns @@ -77,6 +77,14 @@ def _bin_particles_by_radius(data: 'SarracenDataFrame', return rbins, bin_edges +def _get_bin_midpoints(bin_edges: np.ndarray) -> np.ndarray: + """ + Calculate the midpoint of bins given their edges. + """ + + return 0.5 * (bin_edges[1:] - bin_edges[:-1]) + bin_edges[:-1] + + def surface_density(data: 'SarracenDataFrame', r_in: float = None, r_out: float = None, @@ -106,17 +114,17 @@ def surface_density(data: 'SarracenDataFrame', Coordinate system to use to calculate the particle radii. Can be either *spherical* or *cylindrical*. Defaults to *cylindrical*. origin : array-like, optional - The x, y and z position around which to compute radii. Defaults to + The x, y and z centre point around which to compute radii. Defaults to [0, 0, 0]. retbins : bool, optional - Whether to return the bin edges or not. Defaults to False. + Whether to return the midpoints of the bins or not. Defaults to False. Returns ------- array A NumPy array of length bins containing the surface density profile. array, optional - The location of the bin edges. Only returned if *retbins=True*. + The midpoint values of each bin. Only returned if *retbins=True*. Raises ------ @@ -143,7 +151,7 @@ def surface_density(data: 'SarracenDataFrame', sigma = data.groupby(rbins).count().iloc[:, 0] * mass if retbins: - return (sigma / areas).to_numpy(), bin_edges + return (sigma / areas).to_numpy(), _get_bin_midpoints(bin_edges) else: return (sigma / areas).to_numpy() @@ -162,7 +170,7 @@ def _calc_angular_momentum(data: 'SarracenDataFrame', rbins: Series The radial bin to which each particle belongs. origin: list - + The x, y and z centre point around which to compute radii. unit_vector: bool Whether to convert the angular momentum to unit vectors. Default is True. @@ -232,17 +240,17 @@ def angular_momentum(data: 'SarracenDataFrame', Coordinate system to use to calculate the particle radii. Can be either *spherical* or *cylindrical*. Defaults to *cylindrical*. origin : array-like, optional - The x, y and z position around which to compute radii. Defaults to + The x, y and z centre point around which to compute radii. Defaults to [0, 0, 0]. retbins : bool, optional - Whether to return the bin edges or not. Defaults to False. + Whether to return the midpoints of the bins or not. Defaults to False. Returns ------- array A NumPy array of length bins containing the angular momentum profile. array, optional - The location of the bin edges. Only returned if *retbins=True*. + The midpoint values of each bin. Only returned if *retbins=True*. Raises ------ @@ -258,7 +266,7 @@ def angular_momentum(data: 'SarracenDataFrame', Lx, Ly, Lz = Lx.to_numpy(), Ly.to_numpy(), Lz.to_numpy() if retbins: - return Lx, Ly, Lz, bin_edges + return Lx, Ly, Lz, _get_bin_midpoints(bin_edges) else: return Lx, Ly, Lz @@ -276,7 +284,7 @@ def _calc_scale_height(data: 'SarracenDataFrame', rbins: Series The radial bin to which each particle belongs. origin : array-like, optional - The x, y and z position around which to compute radii. Defaults to + The x, y and z centre point around which to compute radii. Defaults to [0, 0, 0]. Returns @@ -301,13 +309,14 @@ def scale_height(data: 'SarracenDataFrame', origin: list = None, retbins: bool = False): """ - Calculates the scale height, H, of the disc. + Calculates the scale height, H/R, of the disc. - The scale height, H, is computed by segmenting the particles into radial + The scale height, H/R, is computed by segmenting the particles into radial bins (rings) and calculating the angular momentum profile of the disc. Each particle takes the dot product of its position vector with the angular momentum vector of its corresponding bin. The standard deviation - of this result per bin yields the scale height profile of the disc. + of this result per bin yields the scale height profile of the disc, which + is divided by the midpoint radius of each bin. Parameters ---------- @@ -324,17 +333,17 @@ def scale_height(data: 'SarracenDataFrame', Coordinate system to use to calculate the particle radii. Can be either *spherical* or *cylindrical*. Defaults to *cylindrical*. origin : array-like, optional - The x, y and z position around which to compute radii. Defaults to + The x, y and z centre point around which to compute radii. Defaults to [0, 0, 0]. retbins : bool, optional - Whether to return the bin edges or not. Defaults to False. + Whether to return the midpoints of the bins or not. Defaults to False. Returns ------- array A NumPy array of length bins scale height, H, profile. array, optional - The location of the bin edges. Only returned if *retbins=True*. + The midpoint values of each bin. Only returned if *retbins=True*. Raises ------ @@ -351,10 +360,11 @@ def scale_height(data: 'SarracenDataFrame', rbins, bin_edges = _bin_particles_by_radius(data, r_in, r_out, bins, geometry, origin) - H = _calc_scale_height(data, rbins, origin).to_numpy() + midpoints = _get_bin_midpoints(bin_edges) + H = _calc_scale_height(data, rbins, origin).to_numpy() / midpoints if retbins: - return H, bin_edges + return H, midpoints else: return H @@ -388,17 +398,17 @@ def honH(data: 'SarracenDataFrame', Coordinate system to use to calculate the particle radii. Can be either *spherical* or *cylindrical*. Defaults to *cylindrical*. origin : array-like, optional - The x, y and z position around which to compute radii. Defaults to + The x, y and z centre point around which to compute radii. Defaults to [0, 0, 0]. retbins : bool, optional - Whether to return the bin edges or not. Defaults to False. + Whether to return the midpoints of the bins or not. Defaults to False. Returns ------- array A NumPy array of length bins containing the /H profile. array, optional - The location of the bin edges. Only returned if *retbins=True*. + The midpoint values of each bin. Only returned if *retbins=True*. Raises ------ @@ -419,6 +429,6 @@ def honH(data: 'SarracenDataFrame', mean_h = data.groupby(rbins)[data.hcol].mean().to_numpy() if retbins: - return mean_h / H, bin_edges + return mean_h / H, _get_bin_midpoints(bin_edges) else: return mean_h / H \ No newline at end of file