From cf296003dd52cb28ff06abc8ec24228b879be363 Mon Sep 17 00:00:00 2001
From: Noe Brucy <noe.brucy@cea.fr>
Date: Thu, 17 Jun 2021 14:50:57 +0200
Subject: [PATCH] [ism] add new file for ism box specific pp

---
 ism.py | 106 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 106 insertions(+)
 create mode 100644 ism.py

diff --git a/ism.py b/ism.py
new file mode 100644
index 0000000..f513474
--- /dev/null
+++ b/ism.py
@@ -0,0 +1,106 @@
+# coding: utf-8
+
+import numpy as np
+import pandas as pd
+from plotter import U
+import postprocessor
+
+mp = 1.4 * 1.66 * 10 ** (-24) * U.g
+z0 = 150 * U.pc
+sink_header = [
+    "Id",
+    "M",
+    "dmf",
+    "x",
+    "y",
+    "z",
+    "vx",
+    "vy",
+    "vz",
+    "rot_period",
+    "lx",
+    "ly",
+    "lz",
+    "acc_rate",
+    "acc_lum",
+    "age",
+    "int_lum",
+    "Teff",
+]
+
+
+def convert_coldens_s(n0):
+    return (np.sqrt(2 * np.pi) * mp * z0 * (n0 * U.cm ** (-3))).express(U.coldens)
+
+
+convert_coldens = np.vectorize(convert_coldens_s)
+
+
+def get_dispersion(dset, i):
+    """
+    Compute dispersion from dset["name"]
+    """
+    vel = dset["vel"][:, i]
+    mass = dset["mass"]
+    mass_tot = np.sum(mass)
+    mean = np.sum(mass * vel) / mass_tot
+    return np.sqrt(np.sum(mass * (vel - mean) ** 2) / mass_tot)
+
+
+def get_sinks(pp):
+    csv_name = f"{pp.path}/output_{pp.num:05}/sink_{pp.num:05}.csv"
+    return pd.read_csv(csv_name, header=None, names=sink_header)
+
+
+def analyze_box(pp):
+    pp.cells["mass"] = postprocessor.mass_func(pp.cells)
+    pp.coldens("z")
+    coldens = pp.get_value("/maps/coldens_z", unit=U.coldens)
+    sinks = get_sinks(pp)
+    sinks["vel"] = np.array([sinks.vx, sinks.vy, sinks.vz])
+    sinks["vel"] *= pp.info["unit_velocity"].express(U.km_s)
+    sinks["mass"] = sinks.M
+    res = {}
+    dirs = ["x", "y", "z"]
+    res["time"] = pp.info["time"] * pp.info["unit_time"].express(U.Myr)
+    for i, dir in enumerate(dirs):
+        res[f"sigma_{dir}"] = get_dispersion(pp.cells, i) * pp.info[
+            "unit_velocity"
+        ].express(U.km_s)
+        res[f"sigma_sinks_{dir}"] = get_dispersion(sinks, i) * pp.info[
+            "unit_velocity"
+        ].express(U.km_s)
+    res["coldens_mean"] = np.mean(coldens)
+    res["n0"] = pp.get_nml("cloud_params/dens0")
+    res["mass"] = np.sum(
+        pp.cells["mass"]
+        * (pp.info["unit_density"] * pp.info["unit_length"] ** 3).express(U.Msun)
+    )
+    res["coldens_initial"] = convert_coldens_s(res["n0"])
+    res["mass_initial"] = res["coldens_initial"] * 1e6
+    res["coldens_mean"] = np.mean(coldens)
+    res["coldens_beam"] = res["mass"] / (pp.info["unit_length"].express(U.pc)) ** 2
+    res["nsink"] = sinks.M.count()
+    res["mass_sink"] = np.sum(sinks.M)
+
+    return res
+
+
+def load_wrapper(pp, fun):
+    """
+    Wrapper to load_unload data and map function
+    """
+    pp.load_cells(keys=["pos", "vel", "dx", "rho"])
+    pp.coldens("z")
+
+    res = fun(pp)
+    pp.unload_cells()
+    pp.unload_parts()
+    del pp.dm
+    del pp.gas
+    del pp.stars
+    return res
+
+
+def allinone(pp):
+    return load_wrapper(pp, analyze_box)