Update galsec from MW

galsec.py

@@ -1,14 +1,21 @@
 # coding: utf-8
 """
     Galaxy kiloparsec extractor
-
+    Noé Brucy 2023
 """
 import numpy as np
 from astropy.table import QTable, hstack
 from astropy import units as u
 from astropy.units.quantity import Quantity
-from scipy.interpolate import griddata
+from collections import defaultdict
-from scipy.fft import fftn
+from numba import jit, prange
+
+
+_atomic_mass = {
+    "H+": 1,
+    "H2": 2,
+    "CO": 28,
+}
 
 
 def vect_r(position: np.array, vector: np.array) -> np.array:
@@ -53,44 +60,92 @@ def vect_phi(position: np.array, vector: np.array):
     return np.einsum("ij,ij -> i", vector[:, :2], uphi)
 
 
+def filter_data(table, bounds):
+
+    mask = np.ones(len(table), dtype=bool)
+    for field in bounds:
+        field_min, field_max = bounds[field]
+        if field_min is not None:
+            mask &= table[field] > field_min
+        if field_max is not None:
+            mask &= table[field] < field_max
+
+    return table[mask]
+
+
+@jit(nopython=True, parallel=True)
 def get_sfr(
-    stars: QTable, time: Quantity[u.Myr], average_on: Quantity[u.Myr] = 30 * u.Myr
+    mass: np.ndarray,
-) -> Quantity[u.Msun / u.year]:
+    birth_time: np.ndarray,
-    """_summary_
+    group_indices: np.ndarray,
+    time: float,
+    average_on: float = 30,
+):
+    """Compute SFR in each group
 
     Parameters
     ----------
-    stars : QTable
+    mass : np.ndarray
-        _description_
+        mass in Msun, dim Ncell, sorted by group
-    time : Quantity[u.Myr],
+    birth_time : np.ndarray
-        time at wich the SFR should be computed
+        birth time in Myr, dim Ncell, sorted by group
-    average_on : Quantity, optional
+    group_indices : np.ndarray
-        compute the sfr on the last `average_on` years, by default 30*u.Myr
+        array of the indices of each bin, dim Ngroup + 1
+    time : float
+        time in Myr
+    average_on : float, optional
+        time sfr is averaged on in Myr, by default 30
 
     Returns
     -------
-    float
+    np.ndarray
-        SFR in Msun/yr
+        sfr for each group in Msun/year, dim Ngroup
     """
-    try:
+
-        mask = (stars["birth_time"] > max(time - average_on, 0 * u.Myr)) & (
+    sfr = np.zeros(len(group_indices) - 1)
-            stars["birth_time"] < time
+    dtime = min(average_on, time)
-        )
+    if dtime > 0:
-        dtime = min(average_on, time).to(u.year)
+        for i in prange(len(group_indices) - 1):
-        if dtime == 0:
+            slice_group = slice(group_indices[i], group_indices[i + 1])
-            sfr = 0.0 * u.Msun / u.year
+            mask = (birth_time[slice_group] > max(time - average_on, 0)) & (
-        else:
+                birth_time[slice_group] < time
-            sfr = np.sum(stars["mass"][mask]) / dtime
+            )
-    except KeyError:
+            sfr[i] = np.sum(mass[slice_group][mask]) / (1e6 * dtime)
-        sfr = 0.0 * u.Msun / u.year
     return sfr
 
 
+@jit(nopython=True, parallel=True)
+def mul_grouped(array: np.ndarray, group_indices: np.ndarray, to_mul: np.ndarray):
+    """Multiply each group by a different value
+
+    Parameters
+    ----------
+    array : np.ndarray
+        array to be multiplied, dim Ncell, sorted by group
+    group_indices : np.ndarray
+        dim Ngroup + 1
+    to_mul : np.ndarray
+        array to multiply with, dim Ngroup
+
+    Returns
+    -------
+    np.ndarray
+        result, dim Ncell
+    """
+    for i in prange(len(group_indices) - 1):
+        slice_group = slice(group_indices[i], group_indices[i + 1])
+        array[slice_group] *= to_mul[i]
+    return array
+
+
 def aggregate(
     grouped_data: QTable,
     weight_field: str = "mass",
-    extensive_fields: str = ["mass", "ek"],
+    extensive_fields: str = [],
     weighted_fields: list = ["velr", "velphi", "velx", "vely", "velz"],
+    compute_dispersions=True,
+    species: list = [],
+    abundance_He=0.1,
 ) -> QTable:
     """Aggregate wisely from grouped data
 
@@ -108,106 +163,138 @@ def aggregate(
     Returns
     -------
     QTable
-        a tabble with the aggregated value for each group
+        a table with the aggregated value for each group
     """
-    assert weight_field in extensive_fields
+    if weight_field not in extensive_fields:
+        extensive_fields.append(weight_field)
+
+    weight_fields_species = []
+    weighted_fields_species = []
+
+    for i, spec in enumerate(species):
+        abundance = grouped_data["chemical_abundances"][:, i]
+        atomic_mass = _atomic_mass[spec]
+        grouped_data[f"{weight_field}_{spec}"] = (
+            grouped_data[weight_field]
+            * abundance
+            * atomic_mass
+            / (1 + 4 * abundance_He)
+        )
+        weight_fields_species.append(f"{weight_field}_{spec}")
 
     for field in weighted_fields:
         # Multiply weighted field by the weight v*m
+        for spec in species:
+            grouped_data[f"{field}_{spec}"] = (
+                grouped_data[field] * grouped_data[f"{weight_field}_{spec}"]
+            )
+            weighted_fields_species.append(f"{field}_{spec}")
         grouped_data[field] *= grouped_data[weight_field]
 
-    # Compute the sum of all fields
+    to_sum = (
-    binned_data = grouped_data[extensive_fields + weighted_fields].groups.aggregate(
+        extensive_fields
-        np.add
+        + weighted_fields
+        + weight_fields_species
+        + weighted_fields_species
     )
+    # Compute the sum of all fields
+    binned_data = grouped_data[to_sum].groups.aggregate(np.add)
 
     for field in weighted_fields:
-        # For weighted field, divided by the total mass
-        # We obtain the weighted mean vmean = 𝚺 (m*v) / 𝚺 m
-        binned_data[field] /= binned_data[weight_field]  
 
-        # We allso compute the weighted dispersion around the weighted mean
+        for spec in ["all"] + species:
-        # Restart from the unweighted data v
+            if spec == "all":
-        grouped_data[field] /= grouped_data[weight_field]
+                weight_field_here = weight_field
-        for i in range(len(grouped_data.groups)):
+                field_here = field
-            # retrieve the indices of each group
+            else:
-            slice = grouped_data.groups.indices[i], grouped_data.groups.indices[i + 1]
+                weight_field_here = f"{weight_field}_{spec}"
-            # Compute the fluctuations wrt the weighted mean (v - vmean)
+                field_here = f"{field}_{spec}"
-            grouped_data[slice[0] : slice[1]][field] -= binned_data[i][field]
 
-        # Compute m * (v - vmean)^2
+            # For weighted field, divided by the total mass
-        grouped_data[field] = grouped_data[weight_field] * grouped_data[field] ** 2
+            # We obtain the weighted mean vmean = 𝚺 (m*v) / 𝚺 m
-        # Compute 𝚺 m * (v - vmean)^2
+            # First we remove 0s to avoid dividing by 0
-        binned_data[f"sigma_{field}"] = grouped_data[field,].groups.aggregate(
+            binned_data[weight_field_here][binned_data[weight_field_here] == 0] = (
-            np.add
+                1.0 * binned_data[weight_field_here].unit
-        )[field]
+            )
-        # Compute sigma = 𝚺 (m * (v - vmean)^2) / 𝚺 m
+            binned_data[field_here] /= binned_data[weight_field_here]
-        binned_data[f"sigma_{field}"] = np.sqrt(
+
-            binned_data[f"sigma_{field}"] / binned_data[weight_field]
+            if compute_dispersions:
-        )
+                # We also compute the weighted dispersion around the weighted mean
+                # First we get m * vmean
+                weight_mean = grouped_data[weight_field_here].value.copy()
+                mul_grouped(
+                    weight_mean,
+                    grouped_data.groups.indices,
+                    binned_data[field_here].value,
+                )
+                # Then we comput (m*v - m*vmean)
+                grouped_data[field_here] -= weight_mean * grouped_data[field_here].unit
+
+                # Compute m * (v - vmean)^2 = (m*v - m*vmean) *  (m*v - m*vmean)  / m
+                # Since we don't want to divide by zero, we first clean them out
+                grouped_data[weight_field_here][
+                    grouped_data[weight_field_here] == 0
+                ] = (1.0 * grouped_data[weight_field_here].unit)
+                grouped_data[field_here] *= (
+                    grouped_data[field_here] / grouped_data[weight_field_here]
+                )
+                if not np.all(np.isfinite(grouped_data[field_here])):
+                    raise FloatingPointError
+                # Compute 𝚺 m * (v - vmean)^2
+                binned_data[f"sigma_{field_here}"] = grouped_data[
+                    field_here,
+                ].groups.aggregate(np.add)[field_here]
+                # Compute sigma = 𝚺 (m * (v - vmean)^2) / 𝚺 m
+                binned_data[f"sigma_{field_here}"] = np.sqrt(
+                    binned_data[f"sigma_{field_here}"] / binned_data[weight_field_here]
+                )
 
     return binned_data
 
 
-def get_bouncing_box_mask(
+def regroup(
-    data: QTable, r: Quantity[u.kpc], phi: Quantity[u.rad], size: Quantity[u.kpc]
+    galex,
+    dataset_key="sectors",
+    bin_key="r",
+    weighted_fields=None,
+    compute_dispersions=True,
 ):
-    x = r * np.cos(phi)
+    dataset = galex.__dict__[dataset_key]
-    y = r * np.sin(phi)
+    result = {}
-    norm_inf = np.maximum(
+    for fluid in galex.fluids:
-        np.abs(data["position"][:, 0] - x), np.abs(data["position"][:, 1] - y)
+        if weighted_fields is None:
-    )
+            weighted_fields_fluid = galex.weighted_fields[fluid]
-    mask = (norm_inf < size / 2) & (np.abs(data["position"][:, 2]) < size / 2)
+        else:
-    return mask
+            weighted_fields_fluid = weighted_fields
 
+        grouped = dataset[fluid].group_by(bin_key)
+        agg_spec = {}
+        for spec in galex.species[fluid]:
+            weighted_fields_spec = [
+                f"{field}_{spec}" for field in weighted_fields_fluid
+            ]
+            agg_spec[spec] = aggregate(
+                grouped,
+                weight_field=f"mass_{spec}",
+                extensive_fields=[],
+                weighted_fields=weighted_fields_spec,
+                compute_dispersions=compute_dispersions,
+            )
+        extensive_fields = []
+        if "sfr" in grouped.keys():
+            extensive_fields.append("sfr")
+        bin_value = grouped[
+            bin_key,
+        ].groups.aggregate(np.fmin)
+        all_values = aggregate(
+            grouped,
+            extensive_fields=extensive_fields,
+            weighted_fields=weighted_fields_fluid,
+            compute_dispersions=compute_dispersions,
+        )
+        result[fluid] = hstack([bin_value, all_values] + list(agg_spec.values()))
+    return result
 
-def regrid(
-    position: Quantity,
-    value: Quantity,
-    resolution: Quantity[u.pc],
-):
-    min_x, max_x = position[:, 0].min(), position[:, 0].max()
-    min_y, max_y = position[:, 1].min(), position[:, 1].max()
-    min_z, max_z = position[:, 2].min(), position[:, 2].max()
-    size = max([max_x - min_x, max_y - min_y, max_z - min_z])
-
-    nb_points = int(np.ceil(((size / resolution).to(u.dimensionless_unscaled).value)))
-    gx, gy, gz = np.mgrid[
-        min_x.value : max_x.value : nb_points * 1j,
-        min_y.value : max_y.value : nb_points * 1j,
-        min_z.value : max_z.value : nb_points * 1j,
-    ]
-
-    grid = griddata(
-        position.value,
-        value.value,
-        (gx, gy, gz),
-        method="nearest",
-    )
-
-    import matplotlib.pyplot as plt
-
-    plt.imshow(
-        grid[:, :, len(grid) // 2].T,
-        origin="lower",
-        extent=[min_x.value, max_x.value, min_y.value, max_y.value],
-    )
-    plt.show()
-
-    return (
-        grid,
-        (gx, gy, gz),
-    )
-
-def fft(
-    position: Quantity,
-    value: Quantity,
-    resolution: Quantity[u.pc],
-    ):
-    
-    grid, (gx, gy, gz) =  regrid(position, value, resolution)
-        
-    fftn(grid, overwrite_x=True)
 
 class Galsec:
     """
@@ -215,22 +302,26 @@ class Galsec:
     """
 
     fluids = ["gas", "stars", "dm"]
+    extensive_fields = defaultdict(
+        lambda: ["mass"],
+        gas=["mass", "volume"],
+    )
+    
+    weighted_fields = defaultdict(lambda: ["velr", "velphi", "velx", "vely", "velz"])
+
     units = {
-        "gas": {
+        "position": u.kpc,
-            "position": u.kpc,
+        "volume": u.kpc**3,
-            "volume": u.kpc**3,
+        "mass": u.Msun,
-            "velocity": u.km / u.s,
+        "density": u.g / u.cm**3,
-            "mass": u.Msun,
+        "velocity": u.km / u.s,
-        },
+        "birth_time": u.Myr,
-        "stars": {
+        "internal_energy": u.km**2 / u.s**2,
-            "position": u.kpc,
-            "velocity": u.km / u.s,
-            "mass": u.Msun,
-            "birth_time": u.Myr,
-        },
-        "dm": {"position": u.kpc, "velocity": u.km / u.s, "mass": u.Msun},
     }
 
+    species = {}
+    abundance_He = 0.1
+
     def __init__(self, data: dict, copy=True) -> None:
         """Initiazise a Galaxasec object
 
@@ -256,9 +347,6 @@ class Galsec:
                 position float array (Ngas, 3) [kpc], centered
                 velocity float array (Ngas, 3) [km/s]
                 mass     float array (Ngas)    [Msun]
-            maps:
-                extent      float list  [xmin, xmax, ymin, ymax]  Coordinates of the edges of the map, centered
-                gas_coldens float array (Nx, Ny) [Msun/pc^2]      Map of column density
         copy : bool, default=True
             Wheter the data should be copied.
         """
@@ -268,11 +356,21 @@ class Galsec:
             self.fluids = data["header"]["fluids"]
 
         for fluid in self.fluids:
-            self.data[fluid] = QTable(data[fluid], units=self.units[fluid], copy=copy)
+            units = {}
+            for key in data[fluid]:
+                if key in self.units:
+                    units[key] = self.units[key]
+            self.data[fluid] = QTable(data[fluid], units=units, copy=copy)
+            self.species[fluid] = []
 
         self.time = data["header"]["time"] * u.Myr
         self.box_size = data["header"]["box_size"] * u.kpc
 
+        if "species" in data["header"]:
+            assert "chemical_abundances" in data["gas"]
+            self.species["gas"] = data["header"]["species"]
+            self.abundance_He = data["header"]["ABHE"]
+
         self.compute_derived_values()
 
     def compute_derived_values(self) -> None:
@@ -283,7 +381,7 @@ class Galsec:
             dset = self.data[fluid]
             dset["r"] = np.sqrt(
                 np.sum(dset["position"][:, :2] ** 2, axis=1)
-            )  # galactic radius
+            )  # galactic radius%run
             dset["phi"] = np.angle(
                 dset["position"][:, 0] + dset["position"][:, 1] * 1j
             )  # galactic longitude
@@ -294,65 +392,80 @@ class Galsec:
                 dset["position"], dset["velocity"]
             )  # azimuthal velocity
             dset["velr"] = vect_r(dset["position"], dset["velocity"])  # radial velocity
-            dset["velx"] = dset["velocity"][:, 0]  # x velocity (TODO this is stupid)
-            dset["vely"] = dset["velocity"][:, 1]  # y velocity
-            dset["velz"] = dset["velocity"][:, 2]  # vertical velocity
-            dset["ek"] = (dset["mass"] * np.sum(dset["velocity"] ** 2, axis=1)).to(
-                u.erg
-            )  # kinetic energy
 
-    def ring_binning(self, delta_r: Quantity[u.kpc]):
+            # Making aliases. No copy is done here.
-        """Add radial bin informations to the dataset
+            dset["x"] = dset["position"][:, 0]
+            dset["y"] = dset["position"][:, 1]
+            dset["z"] = dset["position"][:, 2]
+            dset["velx"] = dset["velocity"][:, 0]
+            dset["vely"] = dset["velocity"][:, 1]
+            dset["velz"] = dset["velocity"][:, 2]
 
-        Parameters
+    def binning(self, bin_deltas):
-        ----------
-        delta_r : Quantity[u.kpc]
-            spacing between two radial bins
-        """
-        for fluid in self.fluids:
-            r_bin = np.trunc(self.data[fluid]["r"] / delta_r)
-            r_bin = (r_bin + 0.5) * delta_r  # Store the middle value
-            self.data[fluid]["r_bin"] = r_bin
-
-    def sector_binning(self, delta_r: Quantity[u.kpc], delta_l: Quantity[u.kpc]):
-        """Add sector bin information to the dataset
-
-        Parameters
-        ----------
-        delta_r : Quantity[u.kpc]
-            spacing between two radial bins
-        delta_l : Quantity[u.kpc]
-            spacing (in spatial units) between two azimuthal bins
-        """
-
-        self.ring_binning(delta_r)
-
-        for fluid in self.fluids:
-            delta_phi = (delta_l / self.data[fluid]["r_bin"]) * u.rad
-            phi_bin = (np.trunc(self.data[fluid]["phi"] / delta_phi) + 0.5) * delta_phi
-            self.data[fluid]["phi_bin"] = phi_bin
-            
-    def cartesian_binning(self, delta_x: Quantity[u.kpc], 
-                          delta_y: Quantity[u.kpc],
-                          delta_z: Quantity[u.kpc]):
-        """Add cartesian bin information to the dataset
-
-        Parameters
-        ----------
-        delta_x : Quantity[u.kpc]
-            spacing between two x bins
-        delta_y : Quantity[u.kpc]
-             spacing between two y bins
-        delta_z : Quantity[u.kpc]
-             spacing between two y bins
-        """
         boxsize = self.box_size
+        keys_binning = []
         for fluid in self.fluids:
-            for i, (delta, name) in enumerate(zip([delta_x, delta_y, delta_z], ["x", "y", "z"])):
+            data = self.data[fluid]  # no copy
-                pos = self.data[fluid]["position"][:, i] + 0.5 * boxsize # stay positive
+            for name in bin_deltas:
-                bin = np.floor(pos / delta)
+                delta = bin_deltas[name]
-                # Store the middle value
+                if delta is not None:
-                self.data[fluid][f"{name}_bin"] = (bin + 0.5) * delta  - 0.5 * boxsize
+                    if name in ["x", "y", "z"]:
+                        # stay positive
+                        pos = data[name] + 0.5 * boxsize
+                        bin = np.floor(pos / delta)
+                        # Store the middle value
+                        data[f"{name}_bin"] = (bin + 0.5) * delta - 0.5 * boxsize
+                    elif name == "r":
+                        r_bin = np.trunc(data["r"] / delta)
+                        r_bin = (r_bin + 0.5) * delta  # Store the middle value
+                        data["r_bin"] = r_bin
+                    elif name == "phi":
+                        delta_phi = (delta / data["r_bin"]) * u.rad
+                        phi_bin = (np.trunc(data["phi"] / delta_phi) + 0.5) * delta_phi
+                        data["phi_bin"] = phi_bin
+                    else:
+                        print(f"Unsupported binning variable {name}")
+                        break
+                    if name not in keys_binning:
+                        keys_binning.append(name)
+        return keys_binning
+
+    def analysis(self, bin_deltas: dict, filter_bounds: dict):
+
+        result = {}
+
+        keys = self.binning(bin_deltas)
+        keys_bin = [key + "_bin" for key in keys]
+
+        for fluid in self.fluids:
+            grouped_data = filter_data(self.data[fluid], filter_bounds).group_by(
+                keys_bin
+            )
+            result[fluid] = hstack(
+                [
+                    grouped_data[keys_bin].groups.aggregate(np.fmin),
+                    aggregate(
+                        grouped_data,
+                        extensive_fields=self.extensive_fields[fluid],
+                        weighted_fields=self.weighted_fields[fluid],
+                        species=self.species[fluid],
+                        abundance_He=self.abundance_He,
+                    ),
+                ]
+            )
+            for key in keys:
+                result[fluid].rename_column(key + "_bin", key)
+
+            if fluid == "stars" and "birth_time" in self.data["stars"].keys():
+                sfr = get_sfr(
+                    grouped_data["mass"].value,
+                    grouped_data["birth_time"].value,
+                    grouped_data.groups.indices,
+                    self.time.value,
+                )
+                result["stars"]["sfr"] = sfr * u.Msun / u.year
+
+        return result
 
     def sector_analysis(
         self,
@@ -363,7 +476,7 @@ class Galsec:
         zmin: Quantity[u.kpc] = -0.5 * u.kpc,
         zmax: Quantity[u.kpc] = 0.5 * u.kpc,
     ):
-        """Compute the aggration of quantities in sectors bins
+        """Compute the aggregation of quantities in sectors bins
 
         Parameters
         ----------
@@ -377,42 +490,10 @@ class Galsec:
             filter out bin beyond that radius, by default 12*u.kpc
         """
 
-        self.sector_binning(delta_r, delta_l)
+        bin_deltas = {"r": delta_r, "phi": delta_l}
-        self.grouped_data = {}
+        filter_bounds = {"r": [rmin, rmax], "z": [zmin, zmax]}
-        self.sectors = {}
+        self.sectors = self.analysis(bin_deltas, filter_bounds)
-
+        return self.sectors
-        for fluid in self.fluids:
-            if fluid == "gas":
-                extensive_fields = ["mass", "ek", "volume"]
-            else:
-                extensive_fields = ["mass", "ek"]
-            filtered_data = self.data[fluid][
-                (self.data[fluid]["r"] > rmin)
-                & (self.data[fluid]["r"] < rmax)
-                & (self.data[fluid]["position"][:, 2] > zmin)
-                & (self.data[fluid]["position"][:, 2] < zmax)
-            ]
-            self.grouped_data[fluid] = filtered_data.group_by(["r_bin", "phi_bin"])
-            self.sectors[fluid] = hstack(
-                [
-                    self.grouped_data[fluid]["r_bin", "phi_bin"].groups.aggregate(
-                        np.fmin
-                    ),
-                    aggregate(
-                        self.grouped_data[fluid], extensive_fields=extensive_fields
-                    ),
-                ]
-            )
-            self.sectors[fluid].rename_column("r_bin", "r")
-            self.sectors[fluid].rename_column("phi_bin", "phi")
-
-        self.sectors["stars"]["sfr"] = (
-            np.zeros(len(self.sectors["stars"]["mass"])) * u.Msun / u.year
-        )
-        for i, group in enumerate(self.grouped_data["stars"].groups):
-            self.sectors["stars"]["sfr"][i] = get_sfr(group, self.time)
-
-        self.sectors["stars"]["sfr"][i] = get_sfr(group, self.time)
 
     def cartesian_analysis(
         self,
@@ -422,7 +503,7 @@ class Galsec:
         zmin: Quantity[u.kpc] = -0.5 * u.kpc,
         zmax: Quantity[u.kpc] = 0.5 * u.kpc,
     ):
-        """Compute the aggration of quantities in cartesian bins
+        """Compute the aggregation of quantities in cartesian bins
 
         Parameters
         ----------
@@ -432,42 +513,10 @@ class Galsec:
              spacing between two y bins
         """
 
-        self.cartesian_binning(delta_x, delta_y, delta_z)
+        bin_deltas = {"x": delta_x, "y": delta_y, "z": delta_z}
-        self.grouped_data = {}
+        filter_bounds = {"z": [zmin, zmax]}
-        self.grid = {}
+        self.grid = self.analysis(bin_deltas, filter_bounds)
-
+        return self.grid
-        for fluid in self.fluids:
-            if fluid == "gas":
-                extensive_fields = ["mass", "ek", "volume"]
-            else:
-                extensive_fields = ["mass", "ek"]
-            filtered_data = self.data[fluid][
-                (self.data[fluid]["position"][:, 2] > zmin)
-                & (self.data[fluid]["position"][:, 2] < zmax)
-            ]
-            self.grouped_data[fluid] = filtered_data.group_by(["x_bin", "y_bin", "z_bin"])
-            self.grid[fluid] = hstack(
-                [
-                    self.grouped_data[fluid]["x_bin", "y_bin", "z_bin"].groups.aggregate(
-                        np.fmin
-                    ),
-                    aggregate(
-                        self.grouped_data[fluid], extensive_fields=extensive_fields
-                    ),
-                ]
-            )
-            self.grid[fluid].rename_column("x_bin", "x")
-            self.grid[fluid].rename_column("y_bin", "y")
-            self.grid[fluid].rename_column("z_bin", "z")
-
-        self.grid["stars"]["sfr"] = (
-            np.zeros(len(self.grid["stars"]["mass"])) * u.Msun / u.year
-        )
-        for i, group in enumerate(self.grouped_data["stars"].groups):
-            self.grid["stars"]["sfr"][i] = get_sfr(group, self.time)
-
-        self.grid["stars"]["sfr"][i] = get_sfr(group, self.time)
-
 
     def ring_analysis(
         self,
@@ -489,34 +538,21 @@ class Galsec:
             filter out bin beyond that radius, by default 30*u.kpc
         """
 
-        self.ring_binning(delta_r)
+        bin_deltas = {"r": delta_r}
-        grouped_data = {}
+        filter_bounds = {"r": [rmin, rmax], "z": [zmin, zmax]}
-        self.rings = {}
+        self.rings = self.analysis(bin_deltas, filter_bounds)
+        return self.rings
 
-        for fluid in self.fluids:
+    def vertical_ring_analysis(
-            if fluid == "gas":
+        self,
-                extensive_fields = ["mass", "ek", "volume"]
+        delta_r: Quantity[u.kpc] = u.kpc,
-            else:
+        delta_z: Quantity[u.kpc] = u.kpc,
-                extensive_fields = ["mass", "ek"]
+        rmin: Quantity[u.kpc] = 1 * u.kpc,
-            filtered_data = self.data[fluid][
+        rmax: Quantity[u.kpc] = 30 * u.kpc,
-                (self.data[fluid]["r"] > rmin)
+        zmin: Quantity[u.kpc] = -0.5 * u.kpc,
-                & (self.data[fluid]["r"] < rmax)
+        zmax: Quantity[u.kpc] = 0.5 * u.kpc,
-                & (self.data[fluid]["position"][:, 2] > zmin)
+    ):
-                & (self.data[fluid]["position"][:, 2] < zmax)
-            ]
-            grouped_data[fluid] = filtered_data.group_by(["r_bin"])
-            self.rings[fluid] = hstack(
-                [
-                    grouped_data[fluid][
-                        "r_bin",
-                    ].groups.aggregate(np.fmin),
-                    aggregate(grouped_data[fluid], extensive_fields=extensive_fields),
-                ]
-            )
-            self.rings[fluid].rename_column("r_bin", "r")
 
-        self.rings["stars"]["sfr"] = (
+        bin_deltas = {"r": delta_r, "z": delta_z}
-            np.zeros(len(self.rings["stars"]["mass"])) * u.Msun / u.year
+        filter_bounds = {"r": [rmin, rmax], "z": [zmin, zmax]}
-        )
+        return self.analysis(bin_deltas, filter_bounds)
-        for i, group in enumerate(grouped_data["stars"].groups):
-            self.rings["stars"]["sfr"][i] = get_sfr(group, self.time)

params_gal.yml

@@ -0,0 +1,96 @@
+plot : # Plot parameters
+    put_title : False                 # Add a title to plot
+
+    # Maps
+    ntick     : 6                     # Number of ticks for maps
+
+    # Overlays
+    vel_red   : 40                    # Take point each vel_red for velocities
+
+    time_fmt  :  "time = {:.3g} {}"   # Time format string, 1st field is time and 2nd is unit
+
+disk: # Disk specific parameters
+    enable    : True                 # Enable specific disk analysis
+    center    : [0.5, 0.5, 0.5]       # Position of the center
+    binning   : "log"                 # Kind of binning (lin : linear, log : logarithmic)
+    mass_star : 1.                    # Mass of the central star
+    nb_bin    : 256                   # Number of bins for averaged quantities
+    bin_in    : 1e-3                  # Outer radius of the inner bin
+    bin_out   : 18                    # Inner radius of the outer bin
+    rmin_pdf  : 1                     # Inner radius for PDF computation
+    rmax_pdf  : 18                    # Outer radius for PDF computation
+
+
+pdf: # parameters for probability density functions
+    nb_bin    : 100                   # Number of bins for the PDF
+    range     : [-1.5, 2.5]           # Range of the PDF (log of fluctuation)
+    xmin_fit  : 0.3                   # Lower boundary of the fit (log of fluctuation)
+    xmax_fit  : 1.5                   # Upper boundary of the fit (log of fluctuation)
+    fit_cut   : 1e-4                  # Exclude value that are < fit_cut * maximum
+
+
+pymses: # Parameters for Pymses reader
+
+    # Source settings
+    variables :  ["rho", "vel", "P"]  # Read these variables
+    part_variables :   ["vel","mass","id","level","epoch"]  # Read these variables
+    order     : '<'                                       # Bit order
+
+
+    # Processing options
+    levelmax : 20                      # Maximal AMR level visited when looking levels
+    fft      : False                   # Quick and dirty rendering using FFT
+    verbose  : True                   # Let pymses write on standart output
+    multiprocessing : True             # Whether to use multiprocessing
+
+    # Camera settings
+    center   :  [0.5, 0.5, 0.5]        # Center of the image
+    zoom     : 4.                      # Zoom  of the image
+    map_size : 2048                    # Size of the computed maps in pixel
+
+    # Filter parameters
+    filter     : False                 # Enable filtering
+    min_coords : [0.35, 0.35, 0.45]
+    max_coords : [0.65, 0.65, 0.55]
+
+input: # Parameters on how to look for input files (= output from Ramses)
+
+    log_prefix     : "run.log"         # Prefix of the log file
+    label_filename : "label.txt"       # Name of the label file
+    nml_filename   : "galaxy.nml"      # name of the namelist file
+    ramses_ism     : False             # If ramses-ism is used
+
+out: # Parameters for post processing
+    tag         : ""                   # Tag for the image
+    interactive : True                 # Interactive mode (keep figures open)
+    save        : True                 # Save the plots on the disk
+    ext         : '.jpeg'              # extension for plots
+    fmt         : ""                   # Format of the output filename for plots
+    # The following keys are accepted
+    # {out} : The output directory (where hdf5  files are also stored)
+    # {run} : Name of the relevant run
+    # {num} : Name of the input file (from Ramses)
+    # {ext} : Extension defined above
+    # {name} : Name of the rule
+    # {tag} : Tag defined above
+    # {nml[nml_key]} : The value of nml_key in the namelist (ex: {nml[amr_params/levelmin]})
+
+
+process: # General setting of the post-processor module
+    verbose : True                     # Give more infos on what is going on
+    num_process : 1                    # Number of forks
+    save_cells : True                  # Save cells structure on disk
+    save_parts : True                  # Save particles on disk
+    unload_cells : True                # Save memory usage
+
+rules: # Specific rules parameters
+    turb_energy_threshold : -1       # Remove invalid data (<0 = no threshold)
+
+
+astrophysix: # Parameters for astrophysix and galactica
+    simu_fmt          : "{tag}_{run}"      # Format of the name of simulation
+    descr_fmt         : "{tag}_{run}"      # Format of the default description
+    # The following keys are accepted
+    # {run} : Name of the relevant run
+    # {tag} : Tag defined above
+    # {nml[nml_key]} : The value of nml_key in the namelist (ex: {nml[amr_params/levelmin]})