From 7429fb8181b50c6d018548e856771d96f37aff19 Mon Sep 17 00:00:00 2001
From: Noe Brucy <noe.brucy@cea.fr>
Date: Wed, 10 May 2023 08:46:13 +0200
Subject: [PATCH] Read boxsize

---
 galsec.py                    | 10 ++++++----
 sectors_extraction_arepo.py  |  4 +++-
 sectors_extraction_ramses.py |  8 +++++---
 3 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/galsec.py b/galsec.py
index eed02f3..9d71eaf 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
+                box_size float [kpc] size of the box
                 fluids   str list list of fluids included
             gas:
                 position float array (Ngas, 3) [kpc], centered
@@ -259,6 +260,7 @@ class Galsec:
             self.data[fluid] = QTable(data[fluid], units=self.units[fluid], copy=copy)
 
         self.time = data["header"]["time"] * u.Myr
+        self.box_size = data["header"]["box_size"] * u.kpc
 
         self.compute_derived_values()
 
@@ -333,13 +335,13 @@ class Galsec:
         delta_z : Quantity[u.kpc]
              spacing between two y bins
         """
+        boxsize = self.box_size
         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]
-                minval = np.floor(np.min(pos) / delta)
-                bin = np.floor((pos - minval) / delta)
+                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  + minval 
+                self.data[fluid][f"{name}_bin"] = (bin + 0.5) * delta  - 0.5 * boxsize
 
     def sector_analysis(
         self,
diff --git a/sectors_extraction_arepo.py b/sectors_extraction_arepo.py
index b58046c..9f78b4c 100644
--- a/sectors_extraction_arepo.py
+++ b/sectors_extraction_arepo.py
@@ -52,7 +52,9 @@ def load_fields(path):
 
     # Create dataset
     data = {
-        "header": {"time": time * UnitTime_in_Myr, "fluids": ["gas", "stars"]},
+        "header": {"time": time * UnitTime_in_Myr, 
+                   "fluids": ["gas", "stars"],
+                   "box_size": box_size * UnitLength_in_kpc },
         "gas": {
             "position": (
                 np.asarray(gas["CenterOfMass"]) - np.array([0.5, 0.5, 0.5]) * box_size
diff --git a/sectors_extraction_ramses.py b/sectors_extraction_ramses.py
index 208c272..f24fea1 100644
--- a/sectors_extraction_ramses.py
+++ b/sectors_extraction_ramses.py
@@ -72,7 +72,9 @@ def load_fields(pp):
 
     # Create dataset
     data = {
-        "header": {"time": pp.time * pp.info["unit_time"].express(U.Myr)},
+        "header": {"time": pp.time * pp.info["unit_time"].express(U.Myr),
+                   "fluids": ["gas", "stars", "dm"],
+                   "box_size": pp.info["unit_length"].express(U.kpc)},
         "gas": {
             "position": cells["position"] * pp.info["unit_length"].express(U.kpc),
             "volume": cells["volume"],
@@ -99,10 +101,10 @@ def load_fields(pp):
     return data
 
 
-if __name__ == "__main__":
+if __name__ == "__main__" and False:
     from snapshotprocessor import SnapshotProcessor
 
     pp = SnapshotProcessor(
-        "/home/nbrucy/ecogal/galturb/F20H_7_1pc_frig", num=80, params="params_gal.yml"
+        "/home/nbrucy/ecogal/galturb/F20H_15_4pc_frig_from_relax", num=50, params="params_gal.yml"
     )
     data = load_fields(pp)