From aaf5c2cb2da2fddef62a27d8ff5ae86499464ad1 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Wed, 10 May 2023 15:37:38 +0200 Subject: [PATCH] add a pipeline for arepo --- galsec.py | 19 +++- ..._extraction_arepo.py => load_data_arepo.py | 10 +-- ...xtraction_ramses.py => load_data_ramses.py | 0 pipeline_MW.py | 90 +++++++++++++++++++ 4 files changed, 107 insertions(+), 12 deletions(-) rename sectors_extraction_arepo.py => load_data_arepo.py (92%) rename sectors_extraction_ramses.py => load_data_ramses.py (100%) create mode 100644 pipeline_MW.py diff --git a/galsec.py b/galsec.py index 9d71eaf..5905ec8 100644 --- a/galsec.py +++ b/galsec.py @@ -113,24 +113,35 @@ def aggregate( assert weight_field in extensive_fields for field in weighted_fields: + # Multiply weighted field by the weight v*m 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 ) for field in weighted_fields: - binned_data[field] /= binned_data[weight_field] # weighted mean + # For weighted field, divided by the total mass + # We obtain the weighted mean vmean = 𝚺 (m*v) / 𝚺 m + binned_data[field] /= binned_data[weight_field] - # Veocity dispersion + # 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] ) @@ -462,7 +473,7 @@ class Galsec: self, delta_r: Quantity[u.kpc] = u.kpc, rmin: Quantity[u.kpc] = 1 * u.kpc, - rmax: Quantity[u.kpc] = 12 * 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, ): @@ -475,7 +486,7 @@ class Galsec: rmin : Quantity[u.kpc], optional filter out bin below that radius, by default 1*u.kpc rmax : Quantity[u.kpc], optional - filter out bin beyond that radius, by default 12*u.kpc + filter out bin beyond that radius, by default 30*u.kpc """ self.ring_binning(delta_r) diff --git a/sectors_extraction_arepo.py b/load_data_arepo.py similarity index 92% rename from sectors_extraction_arepo.py rename to load_data_arepo.py index 9f78b4c..e91e844 100644 --- a/sectors_extraction_arepo.py +++ b/load_data_arepo.py @@ -5,7 +5,7 @@ import h5py from astropy import units as u -def load_fields(path): +def load_fields_arepo(path): """ Parameters ---------- @@ -75,12 +75,6 @@ def load_fields(path): "birth_time": np.asarray(stars["StellarFormationTime"]) * UnitTime_in_Myr, }, } + snap.close() 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/load_data_ramses.py similarity index 100% rename from sectors_extraction_ramses.py rename to load_data_ramses.py diff --git a/pipeline_MW.py b/pipeline_MW.py new file mode 100644 index 0000000..cfd0d58 --- /dev/null +++ b/pipeline_MW.py @@ -0,0 +1,90 @@ +import numpy as np +import matplotlib.pyplot as plt +import h5py +import astropy.units as u, astropy.constants as c + +from load_data_arepo import load_fields_arepo +from galsec import Galsec + +def sfr_history(data, tmin=None, tmax=None, average_time=1, **kwargs): + """History of SFR. Time in Myr. + SFR in Msun/year + Parameters + ---------- + data : _type_ + _description_ + tmin : _type_ + _description_ + tmax : _type_ + _description_ + average_time : _type_ + _description_ + """ + plt.figure(constrained_layout=True, figsize=(5,4)) + if tmin is None: + tmin = np.floor(np.min(data["stars"]["birth_time"])) + if tmax is None: + tmax = np.ceil(np.max(data["stars"]["birth_time"])) + bins = int(np.ceil((tmax - tmin) / average_time)) + tmax = tmin + bins * average_time + plt.hist(data["stars"]["birth_time"], + weights=data["stars"]["mass"] / (average_time*1e6), + bins=bins, range=[tmin, tmax], + histtype="step", **kwargs) + plt.ylabel("SFR [M$_\odot$/yr]") + plt.xlabel("Time [Myr]") + plt.savefig("sfr_history.png") + +def ring_stuff(galsec, delta_r=1): + galsec.ring_analysis(delta_r * u.kpc, rmax=12 * u.kpc) + r_stars = galsec.rings['stars']["r"] + r_gas = galsec.rings['gas']["r"] + surface = np.pi * (delta_r * u.kpc) * r_gas + + plt.figure(constrained_layout=True, figsize=(5,4)) + sfr = galsec.rings['stars']["sfr"] + ssfr = sfr / surface + plt.plot(r_stars, ssfr) + plt.ylabel("SSFR [M$_\odot$/yr/kpc$^2$]") + plt.xlabel("R [kpc]") + plt.savefig("sfr_profile.png") + + + plt.figure(constrained_layout=True, figsize=(5,4)) + velphi = - galsec.rings['gas']["velphi"] + plt.plot(r_gas, velphi) + plt.ylabel(r"$v_\varphi$ [km/s]") + plt.xlabel("R [kpc]") + plt.savefig("rotation_curve.png") + + + plt.figure(constrained_layout=True, figsize=(5,4)) + sigma_velphi = galsec.rings['gas']["sigma_velphi"] + sigma_velr = galsec.rings['gas']["sigma_velr"] + sigma_velz = galsec.rings['gas']["sigma_velz"] + plt.plot(r_gas, sigma_velphi, label=r"$\sigma_{v,\varphi}$") + plt.plot(r_gas, sigma_velr, label=r"$\sigma_{v,r}$") + plt.plot(r_gas, sigma_velz, label=r"$\sigma_{v,z}$") + plt.ylabel(r"$\sigma$ [km/s]") + plt.xlabel("R [kpc]") + plt.legend() + plt.savefig("dispersion_profile.png") + + + +def run_pipeline(path): + + # Galsec analysis + data = load_fields_arepo(path) + galsec = Galsec(data) + + sfr_history(data) + ring_stuff(galsec) + + del data + del galsec + + + +if __name__ == "__main__": + run_pipeline("/home/nbrucy/Travail/Postdoc/Ecogal/MW/junia2_0.25kpc/OUTPUT_SN/snap_400.hdf5") \ No newline at end of file