From e1375ce3367bb0c20f06a62f74b0d05619c7bc96 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Tue, 27 Feb 2024 17:24:00 +0100 Subject: [PATCH] Update galsec from MW --- galsec.py | 606 ++++++++++++++++++++++++++----------------------- params_gal.yml | 96 ++++++++ 2 files changed, 417 insertions(+), 285 deletions(-) create mode 100644 params_gal.yml diff --git a/galsec.py b/galsec.py index 5905ec8..66d537c 100644 --- a/galsec.py +++ b/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 scipy.fft import fftn +from collections import defaultdict +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 -) -> Quantity[u.Msun / u.year]: - """_summary_ + mass: np.ndarray, + birth_time: np.ndarray, + group_indices: np.ndarray, + time: float, + average_on: float = 30, +): + """Compute SFR in each group Parameters ---------- - stars : QTable - _description_ - time : Quantity[u.Myr], - time at wich the SFR should be computed - average_on : Quantity, optional - compute the sfr on the last `average_on` years, by default 30*u.Myr + mass : np.ndarray + mass in Msun, dim Ncell, sorted by group + birth_time : np.ndarray + birth time in Myr, dim Ncell, sorted by group + group_indices : np.ndarray + 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 - SFR in Msun/yr + np.ndarray + sfr for each group in Msun/year, dim Ngroup """ - try: - mask = (stars["birth_time"] > max(time - average_on, 0 * u.Myr)) & ( - stars["birth_time"] < time - ) - dtime = min(average_on, time).to(u.year) - if dtime == 0: - sfr = 0.0 * u.Msun / u.year - else: - sfr = np.sum(stars["mass"][mask]) / dtime - except KeyError: - sfr = 0.0 * u.Msun / u.year + + sfr = np.zeros(len(group_indices) - 1) + dtime = min(average_on, time) + if dtime > 0: + for i in prange(len(group_indices) - 1): + slice_group = slice(group_indices[i], group_indices[i + 1]) + mask = (birth_time[slice_group] > max(time - average_on, 0)) & ( + birth_time[slice_group] < time + ) + sfr[i] = np.sum(mass[slice_group][mask]) / (1e6 * dtime) 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 - binned_data = grouped_data[extensive_fields + weighted_fields].groups.aggregate( - np.add + to_sum = ( + extensive_fields + + 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 - # Restart from the unweighted data v - grouped_data[field] /= grouped_data[weight_field] - for i in range(len(grouped_data.groups)): - # retrieve the indices of each group - slice = grouped_data.groups.indices[i], grouped_data.groups.indices[i + 1] - # Compute the fluctuations wrt the weighted mean (v - vmean) - grouped_data[slice[0] : slice[1]][field] -= binned_data[i][field] - - # Compute m * (v - vmean)^2 - grouped_data[field] = grouped_data[weight_field] * grouped_data[field] ** 2 - # Compute 𝚺 m * (v - vmean)^2 - binned_data[f"sigma_{field}"] = grouped_data[field,].groups.aggregate( - np.add - )[field] - # Compute sigma = 𝚺 (m * (v - vmean)^2) / 𝚺 m - binned_data[f"sigma_{field}"] = np.sqrt( - binned_data[f"sigma_{field}"] / binned_data[weight_field] - ) + for spec in ["all"] + species: + if spec == "all": + weight_field_here = weight_field + field_here = field + else: + weight_field_here = f"{weight_field}_{spec}" + field_here = f"{field}_{spec}" + + # For weighted field, divided by the total mass + # We obtain the weighted mean vmean = 𝚺 (m*v) / 𝚺 m + # First we remove 0s to avoid dividing by 0 + binned_data[weight_field_here][binned_data[weight_field_here] == 0] = ( + 1.0 * binned_data[weight_field_here].unit + ) + binned_data[field_here] /= binned_data[weight_field_here] + + 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( - data: QTable, r: Quantity[u.kpc], phi: Quantity[u.rad], size: Quantity[u.kpc] +def regroup( + galex, + dataset_key="sectors", + bin_key="r", + weighted_fields=None, + compute_dispersions=True, ): - x = r * np.cos(phi) - y = r * np.sin(phi) - norm_inf = np.maximum( - np.abs(data["position"][:, 0] - x), np.abs(data["position"][:, 1] - y) - ) - mask = (norm_inf < size / 2) & (np.abs(data["position"][:, 2]) < size / 2) - return mask + dataset = galex.__dict__[dataset_key] + result = {} + for fluid in galex.fluids: + if weighted_fields is None: + weighted_fields_fluid = galex.weighted_fields[fluid] + else: + 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, - "volume": u.kpc**3, - "velocity": u.km / u.s, - "mass": u.Msun, - }, - "stars": { - "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}, + "position": u.kpc, + "volume": u.kpc**3, + "mass": u.Msun, + "density": u.g / u.cm**3, + "velocity": u.km / u.s, + "birth_time": u.Myr, + "internal_energy": u.km**2 / u.s**2, } + 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]): - """Add radial bin informations to the dataset + # Making aliases. No copy is done here. + 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 - ---------- - 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 - """ + def binning(self, bin_deltas): 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"])): - pos = self.data[fluid]["position"][:, i] + 0.5 * boxsize # stay positive - bin = np.floor(pos / delta) - # Store the middle value - self.data[fluid][f"{name}_bin"] = (bin + 0.5) * delta - 0.5 * boxsize + data = self.data[fluid] # no copy + for name in bin_deltas: + delta = bin_deltas[name] + if delta is not None: + 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,43 +490,11 @@ class Galsec: filter out bin beyond that radius, by default 12*u.kpc """ - self.sector_binning(delta_r, delta_l) - self.grouped_data = {} - self.sectors = {} + bin_deltas = {"r": delta_r, "phi": delta_l} + filter_bounds = {"r": [rmin, rmax], "z": [zmin, zmax]} + 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, delta_x: Quantity[u.kpc] = u.kpc, @@ -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) - self.grouped_data = {} - 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) - + bin_deltas = {"x": delta_x, "y": delta_y, "z": delta_z} + filter_bounds = {"z": [zmin, zmax]} + self.grid = self.analysis(bin_deltas, filter_bounds) + return self.grid 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) - grouped_data = {} - self.rings = {} + bin_deltas = {"r": delta_r} + filter_bounds = {"r": [rmin, rmax], "z": [zmin, zmax]} + self.rings = self.analysis(bin_deltas, filter_bounds) + return self.rings - 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) - ] - 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") + def vertical_ring_analysis( + self, + delta_r: Quantity[u.kpc] = u.kpc, + delta_z: Quantity[u.kpc] = u.kpc, + rmin: Quantity[u.kpc] = 1 * u.kpc, + rmax: Quantity[u.kpc] = 30 * u.kpc, + zmin: Quantity[u.kpc] = -0.5 * u.kpc, + zmax: Quantity[u.kpc] = 0.5 * u.kpc, + ): - self.rings["stars"]["sfr"] = ( - np.zeros(len(self.rings["stars"]["mass"])) * u.Msun / u.year - ) - for i, group in enumerate(grouped_data["stars"].groups): - self.rings["stars"]["sfr"][i] = get_sfr(group, self.time) + bin_deltas = {"r": delta_r, "z": delta_z} + filter_bounds = {"r": [rmin, rmax], "z": [zmin, zmax]} + return self.analysis(bin_deltas, filter_bounds) diff --git a/params_gal.yml b/params_gal.yml new file mode 100644 index 0000000..8bb03ca --- /dev/null +++ b/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]})