diff --git a/galsec.py b/galsec.py index c3db466..a9d83a2 100644 --- a/galsec.py +++ b/galsec.py @@ -229,6 +229,7 @@ class Galsec: A dataset of cells in the following format: header: time float [Myr] time since the start of the simulation + fluids str list list of fluids included gas: position float array (Ngas, 3) [kpc], centered volume float array (Ngas) [pc^3] @@ -251,6 +252,8 @@ class Galsec: """ self.data = {} + if "fluids" in data["header"]: + self.fluids = data["header"]["fluids"] for fluid in self.fluids: self.data[fluid] = QTable(data[fluid], units=self.units[fluid], copy=copy) @@ -406,8 +409,8 @@ class Galsec: filtered_data = self.data[fluid][ (self.data[fluid]["r"] > rmin) & (self.data[fluid]["r"] < rmax) - & (self.data[fluid]["z"] > zmin) - & (self.data[fluid]["z"] < zmax) + & (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( diff --git a/sectors_extraction_arepo.py b/sectors_extraction_arepo.py new file mode 100644 index 0000000..b58046c --- /dev/null +++ b/sectors_extraction_arepo.py @@ -0,0 +1,84 @@ +# coding: utf-8 + +import numpy as np +import h5py +from astropy import units as u + + +def load_fields(path): + """ + Parameters + ---------- + path : _type_ + _description_ + output : _type_ + _description_ + + Returns + ------- + dict + A dataset of cells in the following format: + gas: + position (Ngas, 3) [kpc], centered + volume (Ngas) [pc^3] + velocity (Ngas, 3) [km/s] + mass (Ngas) [Msun] + stars: + position (Nstar, 3) [kpc], centered + velocity (Nstar, 3) [km/s] + mass (Nstar) [Msun] + birth_time (Nstar) [Myr] + dm: + position (Ngas, 3) [kpc], centered + velocity (Ngas, 3) [km/s] + mass (Ngas) [Msun] + maps: + extent (xmin, xmax, ymin, ymax) Coordinates of the edges of the map, centered + gas_coldens (Nx, Ny) [Msun/pc^2], map of column density + """ + + snap = h5py.File(path) + info = snap["Header"].attrs + gas = snap["PartType0"] + stars = snap["PartType4"] + time = info["Time"] + box_size = info["BoxSize"] + UnitLength_in_kpc = info["UnitLength_in_cm"] * u.cm.to(u.kpc) + UnitVelocity_in_km_per_s = info["UnitVelocity_in_cm_per_s"] * u.cm.to(u.km) + UnitTime_in_Myr = ( + info["UnitLength_in_cm"] / info["UnitVelocity_in_cm_per_s"] + ) * u.s.to(u.Myr) + UnitMass_in_Msun = info["UnitMass_in_g"] * u.g.to(u.Msun) + + # Create dataset + data = { + "header": {"time": time * UnitTime_in_Myr, "fluids": ["gas", "stars"]}, + "gas": { + "position": ( + np.asarray(gas["CenterOfMass"]) - np.array([0.5, 0.5, 0.5]) * box_size + ) + * UnitLength_in_kpc, + "volume": (np.asarray(gas["Masses"]) / np.asarray(gas["Density"])) + * UnitLength_in_kpc**3, + "velocity": np.asarray(gas["Velocities"]) * UnitVelocity_in_km_per_s, + "mass": np.asarray(gas["Masses"]) * UnitMass_in_Msun, + }, + "stars": { + "position": ( + np.asarray(stars["Coordinates"]) - np.array([0.5, 0.5, 0.5]) * box_size + ) + * UnitLength_in_kpc, + "velocity": np.asarray(stars["Velocities"]) * UnitVelocity_in_km_per_s, + "mass": np.asarray(stars["Masses"]) * UnitMass_in_Msun, + "birth_time": np.asarray(stars["StellarFormationTime"]) * UnitTime_in_Myr, + }, + } + return data + + +if __name__ == "__main__": + from snapshotprocessor import SnapshotProcessor + + data = load_fields( + "/home/nbrucy/Travail/Postdoc/Ecogal/MW/junia2_0.25kpc/OUTPUT_SN/snap_150.hdf5" + ) diff --git a/sectors_extraction_ramses.py b/sectors_extraction_ramses.py new file mode 100644 index 0000000..208c272 --- /dev/null +++ b/sectors_extraction_ramses.py @@ -0,0 +1,108 @@ +# coding: utf-8 + +import numpy as np +from snapshotprocessor import U +from tables import NoSuchNodeError + + +def load_fields(pp): + """ + Parameters + ---------- + path : _type_ + _description_ + output : _type_ + _description_ + + Returns + ------- + dict + A dataset of cells in the following format: + gas: + position (Ngas, 3) [kpc], centered + volume (Ngas) [pc^3] + velocity (Ngas, 3) [km/s] + mass (Ngas) [Msun] + stars: + position (Nstar, 3) [kpc], centered + velocity (Nstar, 3) [km/s] + mass (Nstar) [Msun] + birth_time (Nstar) [Myr] + dm: + position (Ngas, 3) [kpc], centered + velocity (Ngas, 3) [km/s] + mass (Ngas) [Msun] + maps: + extent (xmin, xmax, ymin, ymax) Coordinates of the edges of the map, centered + gas_coldens (Nx, Ny) [Msun/pc^2], map of column density + """ + + # Load arrays + pp.load_cells(keys=["pos", "vel", "dx", "rho"]) + + try: + pp.load_parts(keys=["pos", "vel", "mass", "epoch"]) + except (KeyError, NoSuchNodeError): + pp.load_parts(keys=["pos", "vel", "mass"]) + + cells = pp.cells + parts = pp.parts + + if "epoch" not in parts: + parts["epoch"] = np.zeros(len(pp.parts["mass"])) + + # Compute extra fields and convert units + for dset in (cells, parts): + dset["position"] = dset["pos"] - np.array([0.5, 0.5, 0.5]) + + cells["mass"] = ( + cells["rho"] + * cells["dx"] ** 3 + * (pp.info["unit_density"] * pp.info["unit_length"] ** 3).express(U.Msun) + ) + + cells["volume"] = cells["dx"] ** 3 * (pp.info["unit_length"] ** 3).express( + U.kpc**3 + ) + + # Separate DM from stars + mass_dm = np.max(parts["mass"]) + mask_dm = parts["mass"] == mass_dm + mask_star = parts["mass"] < mass_dm + + # Create dataset + data = { + "header": {"time": pp.time * pp.info["unit_time"].express(U.Myr)}, + "gas": { + "position": cells["position"] * pp.info["unit_length"].express(U.kpc), + "volume": cells["volume"], + "velocity": cells["vel"] * pp.info["unit_velocity"].express(U.km_s), + "mass": cells["mass"], + }, + "stars": { + "position": parts["position"][mask_star] + * pp.info["unit_length"].express(U.kpc), + "velocity": parts["vel"][mask_star] + * pp.info["unit_velocity"].express(U.km_s), + "mass": parts["mass"][mask_star] * pp.info["unit_mass"].express(U.Msun), + "birth_time": parts["epoch"][mask_star] + * pp.info["unit_time"].express(U.Myr), + }, + "dm": { + "position": parts["position"][mask_dm] + * pp.info["unit_length"].express(U.kpc), + "velocity": parts["vel"][mask_dm] + * pp.info["unit_velocity"].express(U.km_s), + "mass": parts["mass"][mask_dm] * pp.info["unit_mass"].express(U.Msun), + }, + } + return data + + +if __name__ == "__main__": + from snapshotprocessor import SnapshotProcessor + + pp = SnapshotProcessor( + "/home/nbrucy/ecogal/galturb/F20H_7_1pc_frig", num=80, params="params_gal.yml" + ) + data = load_fields(pp)