add a pipeline for arepo

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)

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"
-    )

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")