From 41c1baec44d1683a06854a04a2343bb5bec1c403 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Fri, 18 Feb 2022 22:14:31 +0100
Subject: [PATCH 01/28] added matmul

---
 src/osyris/core/array.py | 23 ++++++++++++++++++++++-
 1 file changed, 22 insertions(+), 1 deletion(-)

diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py
index 1cda9e89..a754a42b 100644
--- a/src/osyris/core/array.py
+++ b/src/osyris/core/array.py
@@ -47,7 +47,10 @@ def __getitem__(self, slice_):
                              unit=self._unit,
                              parent=self._parent,
                              name=self._name)
-        if isinstance(slice_, int) and self.ndim > 1:
+        if self.shape[0] == 1:
+            return out
+        if (isinstance(slice_, int)
+                or isinstance(slice_, np.integer)) and self.ndim > 1:
             return out.reshape(1, len(out))
         return out
 
@@ -253,6 +256,24 @@ def __mul__(self, other):
                                   unit=1.0 * result.units)
         return self.__class__(values=self._array * other, unit=self._unit)
 
+    def __matmul__(self, other):
+        if isinstance(other, self.__class__):
+            scale_l = self._unit.to_base_units()
+            scale_r = other._unit.to_base_units()
+            result = scale_l * scale_r
+            lhs = self._array
+            rhs = other._array * result.magnitude
+            lhs, rhs = self._broadcast(lhs, rhs)
+            return self.__class__(values=np.transpose(lhs @ rhs),
+                                  unit=1.0 * result.units)
+        if isinstance(other, Quantity):
+            scale_l = self._unit.to_base_units()
+            scale_r = other.to_base_units()
+            result = scale_l * scale_r
+            return self.__class__(values=self._array * result.magnitude,
+                                  unit=1.0 * result.units)
+        return self.__class__(values=self._array @ other, unit=self._unit)
+
     def __imul__(self, other):
         if isinstance(other, self.__class__):
             scale_l = self._unit.to_base_units()

From bf66a73dc4816f400f3e1469336484347b8210f2 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Fri, 18 Feb 2022 22:23:17 +0100
Subject: [PATCH 02/28] spatial PR commit

---
 src/osyris/__init__.py         |  1 +
 src/osyris/spatial/__init__.py |  6 +++
 src/osyris/spatial/spatial.py  | 89 ++++++++++++++++++++++++++++++++++
 3 files changed, 96 insertions(+)
 create mode 100644 src/osyris/spatial/__init__.py
 create mode 100644 src/osyris/spatial/spatial.py

diff --git a/src/osyris/__init__.py b/src/osyris/__init__.py
index 0e9f23fd..57e66ecc 100644
--- a/src/osyris/__init__.py
+++ b/src/osyris/__init__.py
@@ -6,3 +6,4 @@
 from .config import config, units
 from .plot import histogram1d, histogram2d, plane, scatter, map, plot
 from .core import Array, Datagroup, Dataset, Plot
+from .spatial import spatial
diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py
new file mode 100644
index 00000000..8db8a97b
--- /dev/null
+++ b/src/osyris/spatial/__init__.py
@@ -0,0 +1,6 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris)
+
+# flake8: noqa
+
+from .spatial import cartesian_to_spherical, spherical_to_cartesian, rotation_matrix
diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py
new file mode 100644
index 00000000..e525c2a4
--- /dev/null
+++ b/src/osyris/spatial/spatial.py
@@ -0,0 +1,89 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris)
+
+import numpy as np
+from ..core import Array
+from ..config import units
+
+
+def rotation_matrix(vec, angle):
+    """
+    Returns 3D rotation matrix of angle 'angle' around rotation vector 'vec'.
+    """
+    if isinstance(vec, list):
+        vec = np.array(vec)
+    vec = vec / np.linalg.norm(vec)
+    R = np.cos(angle) * np.identity(3) + (np.sin(angle)) * np.cross(
+        vec,
+        np.identity(vec.shape[0]) * -1) + (1 - np.cos(angle)) * (np.outer(vec, vec))
+    return R
+
+
+def spherical_to_cartesian(r, theta, phi):
+    """
+    Converts spherical components radius, colatitude and azimuth to x,y,z
+    Returns an osyris array of shape (len(r), 3) containing a stacked xyz
+    """
+    x, y, z = (r * np.cos(phi) * np.sin(theta), r * np.sin(phi) * np.sin(theta),
+               r * np.cos(theta))
+    xyz = np.vstack([x, y, z])
+    return Array(values=np.transpose(xyz).values, unit=r.unit, name="position")
+
+
+def cartesian_to_spherical(position, origin=[0, 0, 0]):
+    """
+    Converts cartesian components in 'position' to spherical components
+    Returns a list of osyris Arrays containing radius, colatitude and azimuth
+    """
+    if isinstance(origin, Array):
+        centered_pos = position - origin
+    else:
+        centered_pos = position - Array(values=origin, unit=centered_pos.unit)
+
+    X = centered_pos[:, 0].values
+    Y = centered_pos[:, 1].values
+    Z = centered_pos[:, 2].values
+
+    radius = np.linalg.norm([X, Y, Z], axis=0)
+    colatitude = np.arctan2(np.linalg.norm([X, Y], axis=0), Z)
+
+    azimuth = np.arctan2(Y, X) + np.pi
+
+    radius = Array(values=radius, unit=position.unit, name='radius')
+    colatitude = Array(values=colatitude, unit=units('rad'), name='colatitude')
+    azimuth = Array(values=azimuth, unit=units('rad'), name='azimuth')
+
+    return radius, colatitude, azimuth
+
+
+def cartesian_to_cylindrical(position, origin=[0, 0, 0]):
+    """
+    Converts cartesian components in 'position' to cylindrical components
+    Returns a list of osyris Arrays containing radius, azimuth and elevation
+    """
+    if isinstance(origin, Array):
+        centered_pos = position - origin
+    else:
+        centered_pos = position - Array(values=origin, unit=centered_pos.unit)
+
+    elevation = centered_pos[:, 2].values
+
+    radius = np.linalg.norm([centered_pos[:, 0].values, centered_pos[:, 1].values],
+                            axis=0)
+    azimuth = np.arctan2(centered_pos[:, 1].values, centered_pos[:, 0].values)
+
+    radius = Array(values=radius, unit=centered_pos.unit, name='radius')
+    azimuth = Array(values=azimuth, unit=units('rad'), name='azimuth')
+    elevation = Array(values=elevation, unit=centered_pos.unit, name='elevation')
+
+    return radius, azimuth, elevation
+
+
+def cylindrical_to_cartesian(radius, azimuth, elevation):
+    """
+    Converts cylindrical components radius, azimuth and elevation to x,y,z
+    Returns an osyris array of shape (len(radius), 3) containing a stacked xyz
+    """
+    x, y, z = (radius * np.cos(azimuth), radius * np.sin(azimuth), elevation)
+    xyz = np.vstack([x, y, z])
+    return Array(values=np.transpose(xyz).values, unit=radius.unit, name="position")

From 7fe084acef6b96068ee005e5b0d61902e0e56f59 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Tue, 15 Mar 2022 16:30:43 +0100
Subject: [PATCH 03/28] array minor change

---
 src/osyris/core/array.py | 121 +++++++++++++++++++++++++++++++++------
 1 file changed, 102 insertions(+), 19 deletions(-)

diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py
index a754a42b..b0b3cb42 100644
--- a/src/osyris/core/array.py
+++ b/src/osyris/core/array.py
@@ -17,27 +17,43 @@ def _comparison_operator(lhs, rhs, op):
     return func(lhs._array, rhs)
 
 
+def _unit_parser(unit):
+    if unit is None:
+        return 1.0 * units.dimensionless
+    elif isinstance(unit, str):
+        return units(unit)
+    elif isinstance(unit, Quantity):
+        return unit
+    elif isinstance(unit, Unit):
+        return 1.0 * unit
+    else:
+        raise TypeError("Unsupported unit type {}".format(type(unit)))
+
+
 class Array:
-    def __init__(self, values=0, unit=None, parent=None, name=""):
+    def __init__(self, values=0, unit=None, parent=None, name="", system="cartesian"):
 
         if isinstance(values, np.ndarray):
             self._array = values
         else:
             self._array = np.asarray(values)
 
-        if unit is None:
-            self._unit = 1.0 * units.dimensionless
-        elif isinstance(unit, str):
-            self._unit = units(unit)
-        elif isinstance(unit, Quantity):
-            self._unit = unit
-        elif isinstance(unit, Unit):
-            self._unit = 1.0 * unit
+        if isinstance(unit, list):
+            units = [_unit_parser(u) for u in unit]
+            self.unit = units
         else:
-            raise TypeError("Unsupported unit type {}".format(type(unit)))
+            self._unit = _unit_parser(unit)
         self._parent = parent
         self._name = name
         self.special_functions = ["sqrt", "power"]
+        coordinate_systems = ["cartesian", "spherical", "cylindrical"]
+        system = system.lower()
+        if system not in coordinate_systems:
+            raise RuntimeError(
+                "Unknown coordinate system keyword '{}'.\nAvailable keywords"
+                " are 'cartesian', 'spherical' and 'cylindrical'.".format(system))
+        else:
+            self._system = system
 
     # def __array__(self):
     #     return self._array
@@ -45,13 +61,20 @@ def __init__(self, values=0, unit=None, parent=None, name=""):
     def __getitem__(self, slice_):
         out = self.__class__(values=self._array[slice_],
                              unit=self._unit,
+                             system=self._system,
                              parent=self._parent,
                              name=self._name)
         if self.shape[0] == 1:
             return out
-        if (isinstance(slice_, int)
-                or isinstance(slice_, np.integer)) and self.ndim > 1:
-            return out.reshape(1, len(out))
+        if isinstance(slice_, (int, np.integer)):
+            if self.ndim > 1:
+                return out.reshape(1, len(out))
+            else:
+                if isinstance(self._unit, list):
+                    out._unit = self._unit[slice_]
+                return out
+        if isinstance(slice_[0], slice) and isinstance(self._unit, list):
+            out._unit = self._unit[slice_[-1]]
         return out
 
     def __len__(self):
@@ -68,9 +91,13 @@ def __str__(self):
             values_str = "Min: " + value_to_string(
                 self.min(use_norm=True).values) + " Max: " + value_to_string(
                     self.max(use_norm=True).values)
-        unit_str = " [{:~}] ".format(self._unit.units)
+        if isinstance(self._unit, list):
+            unit_str = " {} ".format(["{:~}".format(u.units) for u in self._unit])
+        else:
+            unit_str = " [{:~}] ".format(self._unit.units)
         shape_str = str(self._array.shape)
-        return name_str + values_str + unit_str + shape_str
+        system_str = "\nSystem: " + self._system
+        return name_str + values_str + unit_str + shape_str + system_str
 
     def __repr__(self):
         return str(self)
@@ -146,27 +173,81 @@ def name(self, name_):
 
     @property
     def x(self):
-        if self.ndim > 1:
+        if (self.ndim > 1) and (self._system == "cartesian"):
             return self.__class__(values=self._array[:, 0],
                                   unit=self._unit,
                                   parent=self._parent,
                                   name=self._name + "_x")
+        elif self._system != "cartesian":
+            raise AttributeError(
+                "Array '{}' does not have a cartesian basis. Its coordinate system is {}."
+                .format(self._name, self._system))
 
     @property
     def y(self):
-        if self.ndim > 1:
+        if (self.ndim > 1) and (self._system == "cartesian"):
             return self.__class__(values=self._array[:, 1],
                                   unit=self._unit,
                                   parent=self._parent,
                                   name=self._name + "_y")
+        elif self._system != "cartesian":
+            raise AttributeError(
+                "Array '{}' does not have a cartesian basis. Its coordinate system is {}."
+                .format(self._name, self._system))
 
     @property
     def z(self):
-        if self.ndim > 2:
+        if (self.ndim > 2) and (self._system == "cartesian"):
             return self.__class__(values=self._array[:, 2],
                                   unit=self._unit,
                                   parent=self._parent,
                                   name=self._name + "_z")
+        elif self._system != "cartesian":
+            raise AttributeError(
+                "Array '{}' does not have a cartesian basis. Its coordinate system is {}."
+                .format(self._name, self._system))
+
+    @property
+    def r(self):
+        if (self.ndim > 1) and (self._system == "spherical"):
+            unit = self._unit if not isinstance(self._unit, list) else self._unit[0]
+            return self.__class__(values=self._array[:, 0],
+                                  unit=unit,
+                                  parent=self._parent,
+                                  system=self._system,
+                                  name=self._name + "_r")
+        elif self._system != "spherical":
+            raise AttributeError(
+                "Array '{}' does not have a spherical basis. Its coordinate system is {}."
+                .format(self._name, self._system))
+
+    @property
+    def theta(self):
+        if (self.ndim > 1) and (self._system == "spherical"):
+            unit = self._unit if not isinstance(self._unit, list) else self._unit[1]
+            return self.__class__(values=self._array[:, 1],
+                                  unit=unit,
+                                  parent=self._parent,
+                                  system=self._system,
+                                  name=self._name + "_theta")
+        elif self._system != "spherical":
+            raise AttributeError(
+                "Array '{}' does not have a spherical basis. Its coordinate system is {}."
+                .format(self._name, self._system))
+
+    @property
+    def phi(self):
+        if (self.ndim > 2) and (self._system == "spherical"):
+            unit = self._unit if not isinstance(self._unit, list) else self._unit[2]
+            return self.__class__(values=self._array[:, 2],
+                                  unit=unit,
+                                  parent=self._parent,
+                                  system=self._system,
+                                  name=self._name + "_phi")
+        elif self._system != "spherical":
+            raise AttributeError(
+                "Array '{}' does not have a spherical basis. Its coordinate system is {}."
+                .format(self._name, self._system))
 
     @property
     def label(self):
@@ -438,4 +519,6 @@ def max(self, use_norm=False):
         return self.__class__(values=out, unit=self._unit)
 
     def reshape(self, *shape):
-        return self.__class__(values=self._array.reshape(*shape), unit=self._unit)
+        return self.__class__(values=self._array.reshape(*shape),
+                              unit=self._unit,
+                              system=self._system)

From 2a8618d78f0c87d5efd465edbf4d3f0500e80c2b Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 28 Mar 2022 19:06:47 +0200
Subject: [PATCH 04/28] minor stuff

---
 src/osyris/core/array.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py
index b0b3cb42..f2fcfbf8 100644
--- a/src/osyris/core/array.py
+++ b/src/osyris/core/array.py
@@ -96,7 +96,7 @@ def __str__(self):
         else:
             unit_str = " [{:~}] ".format(self._unit.units)
         shape_str = str(self._array.shape)
-        system_str = "\nSystem: " + self._system
+        system_str = " System: " + self._system
         return name_str + values_str + unit_str + shape_str + system_str
 
     def __repr__(self):
@@ -180,7 +180,7 @@ def x(self):
                                   name=self._name + "_x")
         elif self._system != "cartesian":
             raise AttributeError(
-                "Array '{}' does not have a cartesian basis. Its coordinate system is {}."
+                "Array '{}' does not have a cartesian basis. Its coordinate system is {} (attributes are x,y,z)."
                 .format(self._name, self._system))
 
     @property
@@ -192,7 +192,7 @@ def y(self):
                                   name=self._name + "_y")
         elif self._system != "cartesian":
             raise AttributeError(
-                "Array '{}' does not have a cartesian basis. Its coordinate system is {}."
+                "Array '{}' does not have a cartesian basis. Its coordinate system is {} (attributes are x,y,z)."
                 .format(self._name, self._system))
 
     @property
@@ -204,7 +204,7 @@ def z(self):
                                   name=self._name + "_z")
         elif self._system != "cartesian":
             raise AttributeError(
-                "Array '{}' does not have a cartesian basis. Its coordinate system is {}."
+                "Array '{}' does not have a cartesian basis. Its coordinate system is {} (attributes are x,y,z)."
                 .format(self._name, self._system))
 
     @property
@@ -218,7 +218,7 @@ def r(self):
                                   name=self._name + "_r")
         elif self._system != "spherical":
             raise AttributeError(
-                "Array '{}' does not have a spherical basis. Its coordinate system is {}."
+                "Array '{}' does not have a spherical basis. Its coordinate system is {} (attributes are r,theta,phi)."
                 .format(self._name, self._system))
 
     @property
@@ -232,7 +232,7 @@ def theta(self):
                                   name=self._name + "_theta")
         elif self._system != "spherical":
             raise AttributeError(
-                "Array '{}' does not have a spherical basis. Its coordinate system is {}."
+                "Array '{}' does not have a spherical basis. Its coordinate system is {} (attributes are r,theta,phi)."
                 .format(self._name, self._system))
 
     @property
@@ -246,7 +246,7 @@ def phi(self):
                                   name=self._name + "_phi")
         elif self._system != "spherical":
             raise AttributeError(
-                "Array '{}' does not have a spherical basis. Its coordinate system is {}."
+                "Array '{}' does not have a spherical basis. Its coordinate system is {} (attributes are r,theta,phi)."
                 .format(self._name, self._system))
 
     @property

From 86c9e213f5feb3805826004666e5f3a51b097a65 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 25 Apr 2022 12:07:11 +0200
Subject: [PATCH 05/28] temp stash

---
 src/osyris/__init__.py   |   1 +
 src/osyris/core/array.py | 112 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 113 insertions(+)

diff --git a/src/osyris/__init__.py b/src/osyris/__init__.py
index c029893d..4a955b39 100644
--- a/src/osyris/__init__.py
+++ b/src/osyris/__init__.py
@@ -6,3 +6,4 @@
 from .config import config, units
 from .plot import histogram1d, histogram2d, plane, scatter, map, plot
 from .core import Array, Datagroup, Dataset, Plot
+from .spatial import spatial
\ No newline at end of file
diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py
index 9c686849..6ec23dd4 100644
--- a/src/osyris/core/array.py
+++ b/src/osyris/core/array.py
@@ -5,6 +5,7 @@
 from pint.unit import Unit
 from .tools import value_to_string, make_label
 from .. import units
+from .. import spatial as sp
 
 
 def _comparison_operator(lhs, rhs, op):
@@ -162,6 +163,117 @@ def z(self):
                                   parent=self._parent,
                                   name=self._name + "_z")
 
+    def r(self, origin=[0, 0, 0]):
+        radius = False
+        radial_vel = False
+        radial_B = False
+
+        # figure out what kind of computation this is
+        if self.name in ("position", ""):
+            position = self.parent["position"]
+            radius = True
+        elif self.name == "velocity":
+            position = self.parent.parent["amr"]["position"]
+            radial_vel = True
+        elif self.name == "B_field":
+            position = self.parent.parent["amr"]["position"]
+            radial_B = True
+
+
+        # parse origin and compute centered position vector
+        if isinstance(origin, Array):
+            centered_pos = (position - origin).values
+        else:
+            centered_pos = position.values - np.asarray(origin)
+        if radius:
+            r_comp = sp.compute_radius(centered_pos)
+        elif radial_vel:
+            ""
+        elif radial_B:
+            ""
+
+        return self.__class__(name=self.name + "_r", values=r_comp, unit=self._unit)
+
+    def theta(self, origin=[0, 0, 0], basis=None):
+        colatitude = False
+        meridional_vel = False
+        meridional_B = False
+
+        # figure out what kind of computation this is
+        if self.name in ("position", ""):
+            position = self.parent["position"]
+            colatitude = True
+        elif self.name == "velocity":
+            position = self.parent.parent["amr"]["position"]
+            meridional_vel = True
+        elif self.name == "B_field":
+            position = self.parent.parent["amr"]["position"]
+            meridional_B = True
+
+
+        # parse origin and compute centered position vector
+        if isinstance(origin, Array):
+            centered_pos = (position - origin).values
+        else:
+            centered_pos = position.values - np.asarray(origin)
+
+        # rotate coordinates to align with new basis if needed
+        if basis is not None:
+            if isinstance(basis, str):
+                if basis.lower() == "top":
+                    ""
+                elif basis.lower() == "side":
+                    ""
+            else:
+                # normalize vector
+                basis = basis/np.linalg.norm(basis)
+            angle = np.arccos(np.dot([0,0,1], basis))
+            vec = np.cross(basis, [0,0,1])
+            R = sp.rotation_matrix(vec, angle)
+            new_pos = R @ np.transpose(self.parent["position"].values)            
+
+        # compute desired spherical component
+        if colatitude:
+            theta_comp = sp.compute_colatitude(centered_pos)
+        elif meridional_vel:
+            ""
+        elif meridional_B:
+            ""
+
+        return self.__class__(name=self.name + "_theta", values=theta_comp, unit="rad")
+
+    def phi(self, origin=[0, 0, 0], basis=None):
+        azimuth = False
+        azimuthal_vel = False
+        azimuthal_B = False
+
+        # figure out what kind of computation this is
+        if self.name in ("position", ""):
+            position = self.parent["position"]
+            azimuth = True
+        elif self.name == "velocity":
+            position = self.parent.parent["amr"]["position"]
+            radial_vel = True
+        elif self.name == "B_field":
+            position = self.parent.parent["amr"]["position"]
+
+
+        # parse origin and compute centered position vector
+        if isinstance(origin, Array):
+            centered_pos = (position - origin).values
+        else:
+            centered_pos = position.values - np.asarray(origin)
+
+        # compute desired spherical component
+        if azimuth:
+            phi_comp = sp.compute_azimuth(centered_pos)
+        elif azimuthal_vel:
+            ""
+        elif azimuthal_B:
+            ""
+
+        return self.__class__(name=self.name + "_phi", values=phi_comp, unit="rad")
+
     @property
     def label(self):
         return make_label(name=self._name, unit=self._unit.units)

From 5a5acb0da4150dfe1b3df497bb33bdf0d0346e48 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Tue, 31 May 2022 13:34:14 +0200
Subject: [PATCH 06/28] beginning merge

---
 src/osyris/core/vector.py | 253 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 253 insertions(+)
 create mode 100644 src/osyris/core/vector.py

diff --git a/src/osyris/core/vector.py b/src/osyris/core/vector.py
new file mode 100644
index 00000000..a37824c9
--- /dev/null
+++ b/src/osyris/core/vector.py
@@ -0,0 +1,253 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
+import numpy as np
+from pint.quantity import Quantity
+from .base import Base
+from .array import Array
+from .tools import value_to_string
+
+
+def _binary_op(op, lhs, rhs):
+    if isinstance(rhs, (int, float, np.ndarray, Quantity)):
+        rhs = Array(values=rhs)
+    if isinstance(rhs, Array):
+        rhs = lhs.__class__(**{c: rhs for c in lhs._xyz.keys()})
+    if lhs.nvec != rhs.nvec:
+        raise ValueError("Operands do not have the same number of components.")
+
+    return lhs.__class__(
+        **{c: getattr(xyz, op)(getattr(rhs, c))
+           for c, xyz in lhs._xyz.items()})
+
+
+class Vector(Base):
+    def __init__(self, x, y=None, z=None, parent=None, name="", unit=None):
+
+        if isinstance(x, Array):
+            if unit is not None:
+                raise ValueError("Can only set unit when creating Vector from values.")
+            unit = x.unit
+            x = x.values
+            y = self._validate_component(y, x.shape, unit)
+            z = self._validate_component(z, x.shape, unit)
+        self.x = Array(values=x, unit=unit)
+        self.y = Array(values=y, unit=unit) if y is not None else None
+        self.z = Array(values=z, unit=unit) if z is not None else None
+
+        self.parent = parent
+        self.name = name
+
+    def _validate_component(self, array, shape, unit):
+        if array is None:
+            return
+        if array.shape != shape:
+            raise ValueError("The shape of the component does not match the "
+                             "shape of the x component")
+        if array.unit != unit:
+            raise ValueError("The unit of the component does not match the "
+                             "unit of the x component")
+        return array.values
+
+    @property
+    def _xyz(self):
+        out = {'x': self.x}
+        if self.y is not None:
+            out['y'] = self.y
+        if self.z is not None:
+            out['z'] = self.z
+        return out
+
+    def __getitem__(self, slice_):
+        return self.__class__(**{c: xyz[slice_]
+                                 for c, xyz in self._xyz.items()},
+                              parent=self.parent,
+                              name=self._name)
+
+    def __len__(self):
+        return len(self.x)
+
+    def __str__(self):
+        name_str = "'" + self._name + "' "
+        xyz = self._xyz
+        comps_str = ", {" + ",".join(x for x in xyz) + "}"
+        if len(self) == 0:
+            values_str = "Value: " + ", ".join(
+                value_to_string(x.values) for x in xyz.values())
+            unit_str = " [{:~}] ".format(self.unit)
+            shape_str = str(self.shape)
+            return name_str + values_str + unit_str + shape_str + comps_str
+        else:
+            return str(self.norm) + comps_str
+
+    def copy(self):
+        return self.__class__(**{c: xyz.copy()
+                                 for c, xyz in self._xyz.items()},
+                              name=str(self._name))
+
+    @property
+    def norm(self):
+        if (self.y is None) and (self.z is None):
+            return self.x
+        out = self.x.values * self.x.values
+        out += self.y.values * self.y.values
+        if self.z is not None:
+            out += self.z.values * self.z.values
+        return Array(values=np.sqrt(out), unit=self.x.unit, name=self.name)
+
+    @property
+    def unit(self):
+        return self.x.unit
+
+    @unit.setter
+    def unit(self, unit_):
+        self.x.unit = unit_
+        if self.y is not None:
+            self.y.unit = unit_
+        if self.z is not None:
+            self.z.unit = unit_
+
+    @property
+    def ndim(self):
+        return self.x.ndim
+
+    @property
+    def nvec(self):
+        if (self.y is None) and (self.z is None):
+            return 1
+        if self.z is None:
+            return 2
+        return 3
+
+    @property
+    def shape(self):
+        return self.x.shape
+
+    @property
+    def dtype(self):
+        return self.x.dtype
+
+    @property
+    def name(self):
+        return self._name
+
+    @name.setter
+    def name(self, name_):
+        self._name = name_
+        for c, xyz in self._xyz.items():
+            xyz.name = self._name + "_" + c
+
+    def __add__(self, other):
+        return _binary_op("__add__", self, other)
+
+    def __iadd__(self, other):
+        return _binary_op("__iadd__", self, other)
+
+    def __sub__(self, other):
+        return _binary_op("__sub__", self, other)
+
+    def __isub__(self, other):
+        return _binary_op("__isub__", self, other)
+
+    def __mul__(self, other):
+        return _binary_op("__mul__", self, other)
+
+    def __imul__(self, other):
+        return _binary_op("__imul__", self, other)
+
+    def __truediv__(self, other):
+        return _binary_op("__truediv__", self, other)
+
+    def __itruediv__(self, other):
+        return _binary_op("__itruediv__", self, other)
+
+    def __rmul__(self, other):
+        return self * other
+
+    def __rtruediv__(self, other):
+        return np.reciprocal(self / other)
+
+    def __radd__(self, other):
+        return self + other
+
+    def __rsub__(self, other):
+        return -(self - other)
+
+    def __pow__(self, number):
+        return self.__class__(**{c: xyz**number for c, xyz in self._xyz.items()})
+
+    def __neg__(self):
+        return self.__class__(**{c: -xyz for c, xyz in self._xyz.items()})
+
+    def __lt__(self, other):
+        return _binary_op("__lt__", self, other)
+
+    def __le__(self, other):
+        return _binary_op("__le__", self, other)
+
+    def __gt__(self, other):
+        return _binary_op("__gt__", self, other)
+
+    def __ge__(self, other):
+        return _binary_op("__ge__", self, other)
+
+    def __eq__(self, other):
+        return _binary_op("__eq__", self, other)
+
+    def __ne__(self, other):
+        return _binary_op("__ne__", self, other)
+
+    def __and__(self, other):
+        return _binary_op("__and__", self, other)
+
+    def __or__(self, other):
+        return _binary_op("__or__", self, other)
+
+    def __xor__(self, other):
+        return _binary_op("__xor__", self, other)
+
+    def __invert__(self):
+        return np.logical_not(self)
+
+    def to(self, unit):
+        return self.__class__(**{c: xyz.to(unit) for c, xyz in self._xyz.items()})
+
+    def _wrap_numpy(self, func, *args, **kwargs):
+        if isinstance(args[0], (tuple, list)):
+            # Case where we have a sequence of vectors, e.g. `concatenate`
+            out = {
+                c: func(tuple(getattr(a, c) for a in args[0]), *args[1:], **kwargs)
+                for c, xyz in args[0][0]._xyz.items()
+            }
+        elif (len(args) > 1 and isinstance(args[1], self.__class__)):
+            # Case of a binary operation, with two vectors, e.g. `dot`
+            out = {
+                c: func(xyz, getattr(args[1], c), *args[2:], **kwargs)
+                for c, xyz in args[0]._xyz.items()
+            }
+        else:
+            out = {c: func(xyz, *args[1:], **kwargs) for c, xyz in args[0]._xyz.items()}
+        return self.__class__(**out)
+
+    def reshape(self, *shape):
+        return self.__class__(
+            **{c: xyz.reshape(*shape)
+               for c, xyz in self._xyz.items()})
+
+    @property
+    def nbytes(self):
+        return np.sum([xyz.nbytes for xyz in self._xyz.values()])
+
+    def dot(self, other):
+        out = np.zeros(self.shape)
+        for (c1, c2) in zip(self._xyz.values(), other._xyz.values()):
+            out += (c1 * c2).values
+        return Array(values=out, unit=self.unit * other.unit)
+
+    def cross(self, other):
+        x = self.y * other.z
+        x -= self.z * other.y
+        y = self.z * other.x
+        y -= self.x * other.z
+        z = self.x * other.y
+        z -= self.y * other.x
+        return self.__class__(x, y, z)

From 28dff7f23daa35454b9f26b6a2253d5babfb8414 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Thu, 2 Jun 2022 21:24:20 +0200
Subject: [PATCH 07/28] subdomain extraction & coord transforms

---
 src/osyris/core/__init__.py  |   1 +
 src/osyris/core/datagroup.py |   4 +-
 src/osyris/core/subdomain.py | 159 +++++++++++++++++++++++++++++++++++
 src/osyris/core/vector.py    |  25 +++++-
 4 files changed, 186 insertions(+), 3 deletions(-)
 create mode 100644 src/osyris/core/subdomain.py

diff --git a/src/osyris/core/__init__.py b/src/osyris/core/__init__.py
index bb26aeca..049b663b 100644
--- a/src/osyris/core/__init__.py
+++ b/src/osyris/core/__init__.py
@@ -8,3 +8,4 @@
 from .datagroup import Datagroup
 from .dataset import Dataset
 from .plot import Plot
+from .subdomain import Subdomain
\ No newline at end of file
diff --git a/src/osyris/core/datagroup.py b/src/osyris/core/datagroup.py
index 421ca597..5a0f882f 100644
--- a/src/osyris/core/datagroup.py
+++ b/src/osyris/core/datagroup.py
@@ -5,10 +5,10 @@
 
 
 class Datagroup:
-    def __init__(self, data=None, parent=None):
+    def __init__(self, data=None, parent=None, name=""):
         self._container = {}
         self._parent = parent
-        self._name = ""
+        self._name = name
         self.shape = None
         if data is not None:
             for key, array in data.items():
diff --git a/src/osyris/core/subdomain.py b/src/osyris/core/subdomain.py
new file mode 100644
index 00000000..e82a8c3d
--- /dev/null
+++ b/src/osyris/core/subdomain.py
@@ -0,0 +1,159 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
+
+import numpy as np
+from .dataset import Dataset
+from .datagroup import Datagroup
+from .array import Array
+from .vector import Vector
+from .. import spatial
+from .. import units
+
+
+class Subdomain(Dataset):
+    def __init__(self,
+                 dataset,
+                 selection,
+                 dr=None,
+                 dx=None,
+                 dy=None,
+                 dz=None,
+                 origin=None,
+                 basis=None,
+                 dr_L=None):
+
+        self._parent = dataset  # store reference to parent dataset
+        self.meta = dataset.meta
+        self.units = dataset.units
+
+        self._sink_particles = False
+        self._particles = False
+
+        if origin is None:
+            self.origin = Vector(x=0, y=0, z=0, unit=dataset.units['length'].units)
+        else:
+            self.origin = origin
+
+        # find subdomain and extract data
+        valid_comp_amr, valid_comp_sinks, valid_comp_particles = self._determine_subdomain(
+            selection, dr, dx, dy, dz)
+        self._extract_data(valid_comp_amr, valid_comp_sinks, valid_comp_particles)
+
+        # translate positions in subdomain with new origin
+        spatial.change_origin(self, self.origin)
+
+        if basis is None:
+            self.basis = [0, 0, 1]  # assume it's the ramses grid
+        elif isinstance(basis, str):
+            if dr_L is None:
+                raise ValueError(
+                    "Please provide the radius size with which to compute angular momentum (dr_L)"
+                )
+            ang_mom = spatial.get_ang_mom(self, dr_L)
+            if basis.lower() == "top":
+                basis = ang_mom / ang_mom.norm
+            elif basis.lower() == "side":
+                perp_v = Vector(1.0,
+                                1.0,
+                                (-1.0 * (ang_mom.x + ang_mom.y) / ang_mom.z).values,
+                                unit=ang_mom.unit)
+                basis = perp_v / perp_v.norm
+            spatial.change_basis(self, basis)
+        else:
+            spatial.change_basis(self, basis)
+
+    def _determine_subdomain(self, selection, dr, dx, dy, dz):
+        # find indices of data within subdomain for amr dependent groups, sinks and particles
+        centered_pos = self._parent["amr"]["position"] - self.origin
+        if "sink" in self._parent:
+            if len(self._parent["sink"]) > 0:
+                self._sink_particles = True
+                centered_sinks_pos = self._parent["sink"]["position"] - self.origin
+        if "part" in self._parent:
+            if len(self._parent["part"]) > 0:
+                self._particles = True
+                centered_particles_pos = self._parent["part"]["position"] - self.origin
+        if isinstance(selection, str):
+            if selection.lower() == "spherical":
+                if dr is None:
+                    raise ValueError("Please specify a valid selection radius")
+                valid_comp_amr = (centered_pos.norm <= dr).values
+                valid_comp_sinks = (centered_sinks_pos.norm <= dr).values
+                valid_comp_particles = (centered_particles_pos.norm <= dr).values
+            elif selection.lower() == "cubic":
+                if sum(v is not None for v in [dx, dy, dz]) == 3:
+                    # find amr indices
+                    valid_comp_x_amr = (centered_pos.x >= -dx * .5) & (centered_pos.x <=
+                                                                       dx * .5)
+                    valid_comp_y_amr = (centered_pos.y >= -dy * .5) & (centered_pos.y <=
+                                                                       dy * .5)
+                    valid_comp_z_amr = (centered_pos.z >= -dz * .5) & (centered_pos.z <=
+                                                                       dz * .5)
+                    valid_comp_amr = (valid_comp_x_amr & valid_comp_y_amr
+                                      & valid_comp_z_amr).values
+                    # find sink & particles indices
+                    if self._sink_particles:
+                        valid_comp_x_sinks = (centered_sinks_pos.x >= -dx * .5) & (
+                            centered_sinks_pos.x <= dx * .5)
+                        valid_comp_y_sinks = (centered_sinks_pos.y >= -dy * .5) & (
+                            centered_sinks_pos.y <= dy * .5)
+                        valid_comp_z_sinks = (centered_sinks_pos.z >= -dz * .5) & (
+                            centered_sinks_pos.z <= dz * .5)
+                        valid_comp_sinks = (valid_comp_x_sinks & valid_comp_y_sinks
+                                            & valid_comp_z_sinks).values
+                        if not np.any(valid_comp_sinks):
+                            self._sink_particles = False
+                    if self._particles:
+                        valid_comp_x_particles = (centered_particles_pos.x >=
+                                                  -dx * .5) & (centered_particles_pos.x
+                                                               <= dx * .5)
+                        valid_comp_y_particles = (centered_particles_pos.y >=
+                                                  -dy * .5) & (centered_particles_pos.y
+                                                               <= dy * .5)
+                        valid_comp_z_particles = (centered_particles_pos.z >=
+                                                  -dz * .5) & (centered_particles_pos.z
+                                                               <= dz * .5)
+                        valid_comp_particles = (valid_comp_x_particles
+                                                & valid_comp_y_particles
+                                                & valid_comp_z_particles).values
+                        if not np.any(valid_comp_particles):
+                            self._particles = False
+                else:
+                    raise ValueError("Please specify a valid dx, dy and dz")
+            else:
+                raise ValueError(
+                    "Unrecognized string '{}', valid criterions are 'spherical' and 'cubic'"
+                    .format(selection.lower()))
+        elif isinstance(selection, np.array):
+            raise NotImplementedError(
+                "Numpy array indexing for subdomain extraction coming soon...")
+        else:
+            raise ValueError(
+                "Please specify a valid selection criterion (either 'spherical' or 'cubic')"
+            )
+
+        if not np.any(valid_comp_amr):
+            raise ValueError("Empty domain")
+
+        return valid_comp_amr, valid_comp_sinks, valid_comp_particles
+
+    def _extract_data(self, valid_comp_amr, valid_comp_sinks, valid_comp_particles):
+        # proceed to data extraction from parent dataset to subdomain
+        self.groups = {}
+        # extract amr dependent groups
+        for group in ["hydro", "amr", "grav"]:
+            self.groups[group] = Datagroup(parent=self, name=group)
+            for element in self._parent[group]:
+                self.groups[group][element] = self._parent[group][element][
+                    valid_comp_amr]
+        # extract sinks & particles
+        if self._sink_particles:
+            self.groups["sink"] = Datagroup(parent=self, name="sink")
+            for element in self._parent["sink"]:
+                self.groups["sink"][element] = self._parent["sink"][element][
+                    valid_comp_sinks]
+        if self._particles:
+            self.groups["part"] = Datagroup(parent=self, name="part")
+            for element in self._parent["part"]:
+                self.groups["part"][element] = self._parent["part"][element][
+                    valid_comp_particles]
diff --git a/src/osyris/core/vector.py b/src/osyris/core/vector.py
index a37824c9..c0dd906b 100644
--- a/src/osyris/core/vector.py
+++ b/src/osyris/core/vector.py
@@ -5,6 +5,7 @@
 from .base import Base
 from .array import Array
 from .tools import value_to_string
+from .. import spatial
 
 
 def _binary_op(op, lhs, rhs):
@@ -209,7 +210,7 @@ def __invert__(self):
         return np.logical_not(self)
 
     def to(self, unit):
-        return self.__class__(**{c: xyz.to(unit) for c, xyz in self._xyz.items()})
+        return self.__class__(**{c: xyz.to(unit) for c, xyz in self._xyz.items()}, name=self.name)
 
     def _wrap_numpy(self, func, *args, **kwargs):
         if isinstance(args[0], (tuple, list)):
@@ -251,3 +252,25 @@ def cross(self, other):
         z = self.x * other.y
         z -= self.y * other.x
         return self.__class__(x, y, z)
+
+    @property
+    def r(self):
+        if self.name == "position":
+            v = self.norm; v.name = "position_r"
+            return v
+        else:
+            return spatial.get_spherical_components(self, comp='radius')
+
+    @property
+    def theta(self):
+        if self.name == "position":
+            return spatial.get_spherical_colatitude(self)
+        else:
+            return spatial.get_spherical_components(self, comp='colatitude')
+
+    @property
+    def phi(self):
+        if self.name == "position":
+            return spatial.get_spherical_azimuth(self)
+        else:
+            return spatial.get_spherical_components(self, comp='azimuth')

From fe892fc3a763163134ac5f8a57c82de4f0e64882 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Thu, 2 Jun 2022 21:29:38 +0200
Subject: [PATCH 08/28] spatial package added

---
 src/osyris/spatial/__init__.py |   6 ++
 src/osyris/spatial/spatial.py  | 147 +++++++++++++++++++++++++++++++++
 2 files changed, 153 insertions(+)
 create mode 100644 src/osyris/spatial/__init__.py
 create mode 100644 src/osyris/spatial/spatial.py

diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py
new file mode 100644
index 00000000..5b6a590a
--- /dev/null
+++ b/src/osyris/spatial/__init__.py
@@ -0,0 +1,6 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris)
+
+# flake8: noqa
+
+from .spatial import *
\ No newline at end of file
diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py
new file mode 100644
index 00000000..96ffe001
--- /dev/null
+++ b/src/osyris/spatial/spatial.py
@@ -0,0 +1,147 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris)
+
+import numpy as np
+
+
+def change_origin(dataset, new_origin):
+    """
+    Translate all positionnal coordinates to new origin
+    """
+    for g in dataset.groups.keys():
+        for element in dataset[g]:
+            if "position" in element.lower():
+                dataset[g][element] -= new_origin
+    dataset.origin = new_origin
+
+
+def rotation_matrix(vec, angle):
+    """
+    Returns 3D rotation matrix of angle 'angle' around rotation vector 'vec'.
+    """
+    if isinstance(vec, list):
+        vec = np.array(vec)
+    vec = vec / np.linalg.norm(vec)
+    R = np.cos(angle) * np.identity(3) + (np.sin(angle)) * np.cross(
+        vec,
+        np.identity(vec.shape[0]) * -1) + (1 - np.cos(angle)) * (np.outer(vec, vec))
+    return R
+
+
+def cartesian_to_spherical_rotation_matrix(colatitude, azimuth):
+    """
+    Used to convert arbitrary cartesian vectors into spherical vectors
+    """
+
+    R = np.array([[
+        np.sin(colatitude) * np.cos(azimuth),
+        np.sin(colatitude) * np.sin(azimuth),
+        np.cos(colatitude)
+    ],
+                  [
+                      np.cos(colatitude) * np.cos(azimuth),
+                      np.cos(colatitude) * np.sin(azimuth), -np.sin(colatitude)
+                  ], [-np.sin(azimuth), np.cos(azimuth), 0]],
+                 dtype=object)
+    return R
+
+
+def get_ang_mom(subdomain, dr_L):
+    """
+    Compute angular momentum vector in sphere of radius dr_L
+    """
+    sphere = (subdomain["amr"]["position"].norm <= dr_L).values
+    pos = subdomain["amr"]["position"][sphere]
+    mv = subdomain["hydro"]["mass"][sphere] * subdomain["hydro"]["velocity"][sphere]
+    L = np.sum(pos.cross(mv))
+    return L
+
+
+def change_basis(subdomain, new_basis):
+    """
+    Rotates all vectors in dataset to align with vector in 'new_basis'
+    """
+    try:
+        old_basis = subdomain.basis
+    except AttributeError:
+        old_basis = [0, 0, 1]  # assume it's the ramses grid
+    if hasattr(new_basis, 'nvec'):  # if it's a vector
+        new_basis = [new_basis.x.values, new_basis.y.values, new_basis.z.values]
+    rot_angle = np.arccos(
+        np.dot(old_basis, new_basis) /
+        (np.linalg.norm(old_basis) * np.linalg.norm(new_basis)))
+    rot_vector = np.cross(new_basis, old_basis)
+    R = rotation_matrix(rot_vector, rot_angle)
+    for g in subdomain.groups.keys():
+        for element in subdomain[g]:
+            if hasattr(subdomain[g][element], 'nvec'):
+                # all of this will be simplified once matmul is integraded into Vector
+                vector = np.array([
+                    subdomain[g][element].x.values, subdomain[g][element].y.values,
+                    subdomain[g][element].z.values
+                ])
+                vector = R @ vector
+                subdomain[g][element] = subdomain[g][element].__class__(
+                    x=vector[0],
+                    y=vector[1],
+                    z=vector[2],
+                    unit=subdomain[g][element].unit)
+    subdomain.basis = new_basis
+
+
+def get_spherical_colatitude(position_vec):
+    """
+    Compute colatitude from cartesian position vector
+    """
+    X = position_vec.x.values
+    Y = position_vec.y.values
+    Z = position_vec.z.values
+
+    theta = np.arctan2(np.linalg.norm([X, Y], axis=0), Z)
+
+    return position_vec.x.__class__(values=theta, unit='rad', name="position_phi")
+
+
+def get_spherical_azimuth(position_vec):
+    """
+    Compute colatitude from cartesian position vector
+    """
+    X = position_vec.x.values
+    Y = position_vec.y.values
+    phi = np.arctan2(Y, X)
+
+    return position_vec.x.__class__(values=phi, unit='rad', name="position_phi")
+
+
+def get_spherical_components(vec, comp):
+    """
+    Get the spherical components of an arbitrary cartesian vector 'vec'
+    """
+    if not hasattr(vec, 'nvec'):
+        raise ValueError("This is not a vector")
+    if vec.parent.name == "hydro":
+        position = vec.parent.parent["amr"]["position"]
+    else:
+        position = vec.parent["position"]
+
+    theta = position.theta
+    phi = position.phi
+
+    if comp.lower() == "radius":
+        vec_r = np.cos(phi) * np.sin(theta) * vec.x + np.sin(phi) * np.sin(
+            theta) * vec.y + np.cos(theta) * vec.z
+        vec_r.name = vec.name + "_r"
+        return vec_r
+    elif comp.lower() == "colatitude":
+        vec_theta = np.cos(theta) * np.cos(phi) * vec.x + np.cos(theta) * np.sin(
+            phi) * vec.y - np.sin(theta) * vec.z
+        vec_theta.name = vec.name + "_theta"
+        return vec_theta
+    elif comp.lower() == "azimuth":
+        vec_phi = -np.sin(phi) * vec.x + np.cos(phi) * vec.y
+        vec_phi.name = vec.name + "_phi"
+        return vec_phi
+    else:
+        raise ValueError(
+            "Urecognized component keyword '{}'. Valid strings are 'radius', 'colatitude' and 'azimuth'"
+        )

From 9711bea3fbc743305857d5109a9237edc5e414ac Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Thu, 2 Jun 2022 21:41:57 +0200
Subject: [PATCH 09/28] minor yapf

---
 src/osyris/core/vector.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/osyris/core/vector.py b/src/osyris/core/vector.py
index c0dd906b..ea135434 100644
--- a/src/osyris/core/vector.py
+++ b/src/osyris/core/vector.py
@@ -210,7 +210,9 @@ def __invert__(self):
         return np.logical_not(self)
 
     def to(self, unit):
-        return self.__class__(**{c: xyz.to(unit) for c, xyz in self._xyz.items()}, name=self.name)
+        return self.__class__(**{c: xyz.to(unit)
+                                 for c, xyz in self._xyz.items()},
+                              name=self.name)
 
     def _wrap_numpy(self, func, *args, **kwargs):
         if isinstance(args[0], (tuple, list)):
@@ -256,7 +258,8 @@ def cross(self, other):
     @property
     def r(self):
         if self.name == "position":
-            v = self.norm; v.name = "position_r"
+            v = self.norm
+            v.name = "position_r"
             return v
         else:
             return spatial.get_spherical_components(self, comp='radius')

From 7ca0cbca4b9d3bb4c79b135a664f6165ea93cbb8 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Thu, 2 Jun 2022 21:52:26 +0200
Subject: [PATCH 10/28] forgot import of subdomain

---
 src/osyris/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/osyris/__init__.py b/src/osyris/__init__.py
index 67353306..6e215cd9 100644
--- a/src/osyris/__init__.py
+++ b/src/osyris/__init__.py
@@ -5,5 +5,5 @@
 
 from .config import config
 from .units import units
-from .core import Array, Datagroup, Dataset, Plot, Vector
+from .core import Array, Datagroup, Dataset, Plot, Vector, Subdomain
 from .plot import histogram1d, histogram2d, plane, scatter, map, plot

From a0071b3c5d05c72cdfc4181d079c3afac08b29e0 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Tue, 7 Jun 2022 15:00:38 +0200
Subject: [PATCH 11/28] small changes

---
 src/osyris/core/subdomain.py  | 9 +++++++--
 src/osyris/spatial/spatial.py | 2 +-
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/osyris/core/subdomain.py b/src/osyris/core/subdomain.py
index e82a8c3d..48b99530 100644
--- a/src/osyris/core/subdomain.py
+++ b/src/osyris/core/subdomain.py
@@ -73,13 +73,18 @@ def _determine_subdomain(self, selection, dr, dx, dy, dz):
             if len(self._parent["part"]) > 0:
                 self._particles = True
                 centered_particles_pos = self._parent["part"]["position"] - self.origin
+
+        valid_comp_amr, valid_comp_sinks, valid_comp_particles = None, None, None
+
         if isinstance(selection, str):
             if selection.lower() == "spherical":
                 if dr is None:
                     raise ValueError("Please specify a valid selection radius")
                 valid_comp_amr = (centered_pos.norm <= dr).values
-                valid_comp_sinks = (centered_sinks_pos.norm <= dr).values
-                valid_comp_particles = (centered_particles_pos.norm <= dr).values
+                if self._sink_particles:
+                    valid_comp_sinks = (centered_sinks_pos.norm <= dr).values
+                if self._particles:
+                    valid_comp_particles = (centered_particles_pos.norm <= dr).values
             elif selection.lower() == "cubic":
                 if sum(v is not None for v in [dx, dy, dz]) == 3:
                     # find amr indices
diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py
index 96ffe001..368539ab 100644
--- a/src/osyris/spatial/spatial.py
+++ b/src/osyris/spatial/spatial.py
@@ -143,5 +143,5 @@ def get_spherical_components(vec, comp):
         return vec_phi
     else:
         raise ValueError(
-            "Urecognized component keyword '{}'. Valid strings are 'radius', 'colatitude' and 'azimuth'"
+            "Urecognized component keyword '{}'. Valid strings are 'radius', 'colatitude' and 'azimuth'".format(comp)
         )

From 8a26e2bf32b5152037663aa16d1a8cf120c37181 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Tue, 7 Jun 2022 15:01:21 +0200
Subject: [PATCH 12/28] cylindrical transforms

---
 src/osyris/core/vector.py     |  9 +++++++++
 src/osyris/spatial/spatial.py | 34 +++++++++++++++++++++++++++++++---
 2 files changed, 40 insertions(+), 3 deletions(-)

diff --git a/src/osyris/core/vector.py b/src/osyris/core/vector.py
index ea135434..f0f46acd 100644
--- a/src/osyris/core/vector.py
+++ b/src/osyris/core/vector.py
@@ -277,3 +277,12 @@ def phi(self):
             return spatial.get_spherical_azimuth(self)
         else:
             return spatial.get_spherical_components(self, comp='azimuth')
+
+    @property
+    def cyl_r(self):
+        if self.name == "position":
+            v = Array(values=np.sqrt(self.x.values**2 + self.y.values**2), unit=self.unit)
+            v.name = "position_cyl_r"
+            return v
+        else:
+            return spatial.get_cylindrical_radial_component(self)
\ No newline at end of file
diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py
index 368539ab..2c699242 100644
--- a/src/osyris/spatial/spatial.py
+++ b/src/osyris/spatial/spatial.py
@@ -46,6 +46,17 @@ def cartesian_to_spherical_rotation_matrix(colatitude, azimuth):
     return R
 
 
+def cartesian_to_cylindrical_rotation_matrix(azimuth):
+    """
+    Used to convert arbitrary cartesian vectors into cylindrical vectors
+    """
+
+    R = np.array([[np.cos(azimuth), np.sin(azimuth), 0],
+                  [-np.sin(azimuth), np.cos(azimuth), 0], [0, 0, 1]],
+                 dtype=object)
+    return R
+
+
 def get_ang_mom(subdomain, dr_L):
     """
     Compute angular momentum vector in sphere of radius dr_L
@@ -104,7 +115,7 @@ def get_spherical_colatitude(position_vec):
 
 def get_spherical_azimuth(position_vec):
     """
-    Compute colatitude from cartesian position vector
+    Compute azimuth from cartesian position vector
     """
     X = position_vec.x.values
     Y = position_vec.y.values
@@ -113,6 +124,23 @@ def get_spherical_azimuth(position_vec):
     return position_vec.x.__class__(values=phi, unit='rad', name="position_phi")
 
 
+def get_cylindrical_radial_component(vec):
+    """
+    Get the radial component of an arbitrary cartesian vector 'vec' in cylindrical coords
+    """
+    if not hasattr(vec, 'nvec'):
+        raise ValueError("This is not a vector")
+    if vec.parent.name == "hydro":
+        position = vec.parent.parent["amr"]["position"]
+    else:
+        position = vec.parent["position"]
+
+    phi = position.phi
+    vec_r = np.cos(phi)*vec.x + np.sin(phi)*vec.y
+    vec_r.name = vec.name + "_cyl_r"
+
+    return vec_r
+
 def get_spherical_components(vec, comp):
     """
     Get the spherical components of an arbitrary cartesian vector 'vec'
@@ -143,5 +171,5 @@ def get_spherical_components(vec, comp):
         return vec_phi
     else:
         raise ValueError(
-            "Urecognized component keyword '{}'. Valid strings are 'radius', 'colatitude' and 'azimuth'".format(comp)
-        )
+            "Urecognized component keyword '{}'. Valid strings are 'radius', 'colatitude' and 'azimuth'"
+            .format(comp))

From a4800cbdf960699e0ee1a548aea2b994239c3f18 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 13 Jun 2022 16:49:00 +0200
Subject: [PATCH 13/28] slicing PR. This is subdomain extraction

---
 docs/basics.ipynb               |   2 +-
 docs/data_structures.ipynb      |   2 +-
 docs/loading_data.ipynb         |   4 +-
 docs/plotting_maps.ipynb        |   2 +-
 docs/recipes.ipynb              |   2 +-
 setup.cfg                       |   2 +-
 src/osyris/__init__.py          |   2 +-
 src/osyris/config/defaults.py   |  48 +++------
 src/osyris/core/__init__.py     |   3 +-
 src/osyris/core/array.py        |   8 ++
 src/osyris/core/dataset.py      |  12 ++-
 src/osyris/core/subdomain.py    | 164 ----------------------------
 src/osyris/io/sink.py           |   2 +-
 src/osyris/spatial/spatial.py   | 186 +++++---------------------------
 src/osyris/units/__init__.py    |   7 ++
 src/osyris/units/library.py     |  45 ++++++++
 src/osyris/{ => units}/units.py |   2 +-
 test/test_array.py              |  38 +++++--
 test/test_vector.py             |  23 ++++
 19 files changed, 171 insertions(+), 383 deletions(-)
 delete mode 100644 src/osyris/core/subdomain.py
 create mode 100644 src/osyris/units/__init__.py
 create mode 100644 src/osyris/units/library.py
 rename src/osyris/{ => units}/units.py (96%)

diff --git a/docs/basics.ipynb b/docs/basics.ipynb
index 37e62e8e..b37c27ef 100644
--- a/docs/basics.ipynb
+++ b/docs/basics.ipynb
@@ -161,7 +161,7 @@
    "outputs": [],
    "source": [
     "ind = np.argmax(data[\"hydro\"][\"density\"])\n",
-    "center = data[\"amr\"][\"position\"][ind.values]\n",
+    "center = data[\"amr\"][\"position\"][ind]\n",
     "center"
    ]
   },
diff --git a/docs/data_structures.ipynb b/docs/data_structures.ipynb
index a3395c90..44605cc9 100644
--- a/docs/data_structures.ipynb
+++ b/docs/data_structures.ipynb
@@ -303,7 +303,7 @@
    "outputs": [],
    "source": [
     "ind = np.argmax(data[\"hydro\"][\"density\"])\n",
-    "center = data[\"amr\"][\"position\"][ind.values]\n",
+    "center = data[\"amr\"][\"position\"][ind]\n",
     "center"
    ]
   },
diff --git a/docs/loading_data.ipynb b/docs/loading_data.ipynb
index 4ba161fd..76c87ec4 100644
--- a/docs/loading_data.ipynb
+++ b/docs/loading_data.ipynb
@@ -119,7 +119,7 @@
    "source": [
     "osyris.map({\"data\": data[\"hydro\"][\"density\"], \"norm\": \"log\"},\n",
     "           dx=1000 * osyris.units('au'),\n",
-    "           origin=data[\"amr\"][\"position\"][np.argmax(data[\"hydro\"][\"density\"]).values],\n",
+    "           origin=data[\"amr\"][\"position\"][np.argmax(data[\"hydro\"][\"density\"])],\n",
     "           direction='z')"
    ]
   },
@@ -252,7 +252,7 @@
    "source": [
     "osyris.map({\"data\": data[\"hydro\"][\"density\"], \"norm\": \"log\"},\n",
     "           dx=2000 * osyris.units(\"au\"),\n",
-    "           origin=data[\"amr\"][\"position\"][np.argmax(data[\"hydro\"][\"density\"]).values],\n",
+    "           origin=data[\"amr\"][\"position\"][np.argmax(data[\"hydro\"][\"density\"])],\n",
     "           direction='z')"
    ]
   },
diff --git a/docs/plotting_maps.ipynb b/docs/plotting_maps.ipynb
index 3fad61b1..fee25e53 100644
--- a/docs/plotting_maps.ipynb
+++ b/docs/plotting_maps.ipynb
@@ -27,7 +27,7 @@
     "path = \"osyrisdata/starformation\"\n",
     "data = osyris.Dataset(8, path=path).load()\n",
     "ind = np.argmax(data[\"hydro\"][\"density\"])\n",
-    "center = data[\"amr\"][\"position\"][ind.values]"
+    "center = data[\"amr\"][\"position\"][ind]"
    ]
   },
   {
diff --git a/docs/recipes.ipynb b/docs/recipes.ipynb
index 2604209b..d0884909 100644
--- a/docs/recipes.ipynb
+++ b/docs/recipes.ipynb
@@ -94,7 +94,7 @@
    "source": [
     "# Find center\n",
     "ind = np.argmax(data8[\"hydro\"][\"density\"])\n",
-    "center = data8[\"amr\"][\"position\"][ind.values]\n",
+    "center = data8[\"amr\"][\"position\"][ind]\n",
     "\n",
     "dx = 2000 * osyris.units(\"au\")\n",
     "# Extract density slices by copying data into structures\n",
diff --git a/setup.cfg b/setup.cfg
index e349711f..3a8dee35 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = osyris
-version = 2.9.0
+version = 2.9.1
 author = Neil Vaytet
 author_email = neil.vaytet@esss.se
 description = A package to visualize AMR data from the RAMSES code
diff --git a/src/osyris/__init__.py b/src/osyris/__init__.py
index 6e215cd9..67353306 100644
--- a/src/osyris/__init__.py
+++ b/src/osyris/__init__.py
@@ -5,5 +5,5 @@
 
 from .config import config
 from .units import units
-from .core import Array, Datagroup, Dataset, Plot, Vector, Subdomain
+from .core import Array, Datagroup, Dataset, Plot, Vector
 from .plot import histogram1d, histogram2d, plane, scatter, map, plot
diff --git a/src/osyris/config/defaults.py b/src/osyris/config/defaults.py
index 3268b4dc..22981703 100644
--- a/src/osyris/config/defaults.py
+++ b/src/osyris/config/defaults.py
@@ -4,7 +4,6 @@
 Define default values so that you don't have to specify them every time.
 """
 
-from collections import defaultdict
 from math import sqrt, pi
 
 
@@ -36,68 +35,45 @@ def configure_units(units, unit_d, unit_l, unit_t):
     temperature = 1.0 * units("K")
     grav_potential = velocity**2
 
-    library = defaultdict(lambda: 1.0 * units('dimensionless'))
-
-    library.update({
+    library = {
         'unit_d': unit_d,
         'unit_l': unit_l,
         'unit_t': unit_t,
         'density': density,
         'velocity': velocity,
-        'velocity_x': velocity,
-        'velocity_y': velocity,
-        'velocity_z': velocity,
+        'velocity_*': velocity,
         'momentum': momentum,
-        'momentum_x': momentum,
-        'momentum_y': momentum,
-        'momentum_z': momentum,
+        'momentum_*': momentum,
         'magnetic_field': magnetic_field,
         'B_left': magnetic_field,
-        'B_left_x': magnetic_field,
-        'B_left_y': magnetic_field,
-        'B_left_z': magnetic_field,
+        'B_left_*': magnetic_field,
         'B_right': magnetic_field,
-        'B_right_x': magnetic_field,
-        'B_right_y': magnetic_field,
-        'B_right_z': magnetic_field,
+        'B_right_*': magnetic_field,
         'B_field': magnetic_field,
-        'B_field_x': magnetic_field,
-        'B_field_y': magnetic_field,
-        'B_field_z': magnetic_field,
-        'B_x_left': magnetic_field,
-        'B_y_left': magnetic_field,
-        'B_z_left': magnetic_field,
-        'B_x_right': magnetic_field,
-        'B_y_right': magnetic_field,
-        'B_z_right': magnetic_field,
+        'B_field_*': magnetic_field,
+        'B_*_left': magnetic_field,
+        'B_*_right': magnetic_field,
         'acceleration': acceleration,
         'grav_acceleration': acceleration,
-        'grav_acceleration_x': acceleration,
-        'grav_acceleration_y': acceleration,
-        'grav_acceleration_z': acceleration,
+        'grav_acceleration_*': acceleration,
         'grav_potential': grav_potential,
         'energy': energy,
         'internal_energy': energy,
         'thermal_pressure': energy,
         'pressure': energy,
         'radiative_energy': energy,
-        'radiative_energy_1': energy,
+        'radiative_energy_*': energy,
         'time': time,
         'length': length,
         'x': length,
         'y': length,
         'z': length,
-        'xyz_x': length,
-        'xyz_y': length,
-        'xyz_z': length,
         'position': length,
-        'position_x': length,
-        'position_y': length,
-        'position_z': length,
+        'position_*': length,
         'dx': length,
         'mass': mass,
         'temperature': temperature
-    })
+    }
     return library
 
 
diff --git a/src/osyris/core/__init__.py b/src/osyris/core/__init__.py
index 049b663b..e5ab547a 100644
--- a/src/osyris/core/__init__.py
+++ b/src/osyris/core/__init__.py
@@ -7,5 +7,4 @@
 from .vector import Vector
 from .datagroup import Datagroup
 from .dataset import Dataset
-from .plot import Plot
-from .subdomain import Subdomain
\ No newline at end of file
+from .plot import Plot
\ No newline at end of file
diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py
index 01edbcc0..23423424 100644
--- a/src/osyris/core/array.py
+++ b/src/osyris/core/array.py
@@ -48,6 +48,14 @@ def __init__(self, values, unit=None, parent=None, name=""):
         self.name = name
 
     def __getitem__(self, slice_):
+        if isinstance(slice_, Base):
+            if not isinstance(slice_, self.__class__):
+                raise ValueError(
+                    "Cannot slice using a Vector, only Array is supported.")
+            if slice_.dtype not in ('int32', 'int64', bool):
+                raise TypeError(
+                    "Dtype of the Array must be integer or bool for slicing.")
+            slice_ = slice_.values
         return self.__class__(values=self._array[slice_],
                               unit=self.unit,
                               parent=self.parent,
diff --git a/src/osyris/core/dataset.py b/src/osyris/core/dataset.py
index e204dd80..e8847375 100644
--- a/src/osyris/core/dataset.py
+++ b/src/osyris/core/dataset.py
@@ -6,7 +6,7 @@
 from ..io import Loader
 from .datagroup import Datagroup
 from .tools import bytes_to_human_readable
-from .. import units
+from ..units import units, UnitsLibrary
 
 
 class Dataset:
@@ -111,7 +111,9 @@ def update(self, d):
             self[key] = value
 
     def set_units(self):
-        self.units = config.configure_units(units=units,
-                                            unit_d=self.meta['unit_d'],
-                                            unit_l=self.meta['unit_l'],
-                                            unit_t=self.meta['unit_t'])
+        self.units = UnitsLibrary(library=config.configure_units(
+            units=units,
+            unit_d=self.meta['unit_d'],
+            unit_l=self.meta['unit_l'],
+            unit_t=self.meta['unit_t']),
+                                  default_unit=1.0 * units(''))
diff --git a/src/osyris/core/subdomain.py b/src/osyris/core/subdomain.py
deleted file mode 100644
index 48b99530..00000000
--- a/src/osyris/core/subdomain.py
+++ /dev/null
@@ -1,164 +0,0 @@
-# SPDX-License-Identifier: BSD-3-Clause
-# Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
-
-import numpy as np
-from .dataset import Dataset
-from .datagroup import Datagroup
-from .array import Array
-from .vector import Vector
-from .. import spatial
-from .. import units
-
-
-class Subdomain(Dataset):
-    def __init__(self,
-                 dataset,
-                 selection,
-                 dr=None,
-                 dx=None,
-                 dy=None,
-                 dz=None,
-                 origin=None,
-                 basis=None,
-                 dr_L=None):
-
-        self._parent = dataset  # store reference to parent dataset
-        self.meta = dataset.meta
-        self.units = dataset.units
-
-        self._sink_particles = False
-        self._particles = False
-
-        if origin is None:
-            self.origin = Vector(x=0, y=0, z=0, unit=dataset.units['length'].units)
-        else:
-            self.origin = origin
-
-        # find subdomain and extract data
-        valid_comp_amr, valid_comp_sinks, valid_comp_particles = self._determine_subdomain(
-            selection, dr, dx, dy, dz)
-        self._extract_data(valid_comp_amr, valid_comp_sinks, valid_comp_particles)
-
-        # translate positions in subdomain with new origin
-        spatial.change_origin(self, self.origin)
-
-        if basis is None:
-            self.basis = [0, 0, 1]  # assume it's the ramses grid
-        elif isinstance(basis, str):
-            if dr_L is None:
-                raise ValueError(
-                    "Please provide the radius size with which to compute angular momentum (dr_L)"
-                )
-            ang_mom = spatial.get_ang_mom(self, dr_L)
-            if basis.lower() == "top":
-                basis = ang_mom / ang_mom.norm
-            elif basis.lower() == "side":
-                perp_v = Vector(1.0,
-                                1.0,
-                                (-1.0 * (ang_mom.x + ang_mom.y) / ang_mom.z).values,
-                                unit=ang_mom.unit)
-                basis = perp_v / perp_v.norm
-            spatial.change_basis(self, basis)
-        else:
-            spatial.change_basis(self, basis)
-
-    def _determine_subdomain(self, selection, dr, dx, dy, dz):
-        # find indices of data within subdomain for amr dependent groups, sinks and particles
-        centered_pos = self._parent["amr"]["position"] - self.origin
-        if "sink" in self._parent:
-            if len(self._parent["sink"]) > 0:
-                self._sink_particles = True
-                centered_sinks_pos = self._parent["sink"]["position"] - self.origin
-        if "part" in self._parent:
-            if len(self._parent["part"]) > 0:
-                self._particles = True
-                centered_particles_pos = self._parent["part"]["position"] - self.origin
-
-        valid_comp_amr, valid_comp_sinks, valid_comp_particles = None, None, None
-
-        if isinstance(selection, str):
-            if selection.lower() == "spherical":
-                if dr is None:
-                    raise ValueError("Please specify a valid selection radius")
-                valid_comp_amr = (centered_pos.norm <= dr).values
-                if self._sink_particles:
-                    valid_comp_sinks = (centered_sinks_pos.norm <= dr).values
-                if self._particles:
-                    valid_comp_particles = (centered_particles_pos.norm <= dr).values
-            elif selection.lower() == "cubic":
-                if sum(v is not None for v in [dx, dy, dz]) == 3:
-                    # find amr indices
-                    valid_comp_x_amr = (centered_pos.x >= -dx * .5) & (centered_pos.x <=
-                                                                       dx * .5)
-                    valid_comp_y_amr = (centered_pos.y >= -dy * .5) & (centered_pos.y <=
-                                                                       dy * .5)
-                    valid_comp_z_amr = (centered_pos.z >= -dz * .5) & (centered_pos.z <=
-                                                                       dz * .5)
-                    valid_comp_amr = (valid_comp_x_amr & valid_comp_y_amr
-                                      & valid_comp_z_amr).values
-                    # find sink & particles indices
-                    if self._sink_particles:
-                        valid_comp_x_sinks = (centered_sinks_pos.x >= -dx * .5) & (
-                            centered_sinks_pos.x <= dx * .5)
-                        valid_comp_y_sinks = (centered_sinks_pos.y >= -dy * .5) & (
-                            centered_sinks_pos.y <= dy * .5)
-                        valid_comp_z_sinks = (centered_sinks_pos.z >= -dz * .5) & (
-                            centered_sinks_pos.z <= dz * .5)
-                        valid_comp_sinks = (valid_comp_x_sinks & valid_comp_y_sinks
-                                            & valid_comp_z_sinks).values
-                        if not np.any(valid_comp_sinks):
-                            self._sink_particles = False
-                    if self._particles:
-                        valid_comp_x_particles = (centered_particles_pos.x >=
-                                                  -dx * .5) & (centered_particles_pos.x
-                                                               <= dx * .5)
-                        valid_comp_y_particles = (centered_particles_pos.y >=
-                                                  -dy * .5) & (centered_particles_pos.y
-                                                               <= dy * .5)
-                        valid_comp_z_particles = (centered_particles_pos.z >=
-                                                  -dz * .5) & (centered_particles_pos.z
-                                                               <= dz * .5)
-                        valid_comp_particles = (valid_comp_x_particles
-                                                & valid_comp_y_particles
-                                                & valid_comp_z_particles).values
-                        if not np.any(valid_comp_particles):
-                            self._particles = False
-                else:
-                    raise ValueError("Please specify a valid dx, dy and dz")
-            else:
-                raise ValueError(
-                    "Unrecognized string '{}', valid criterions are 'spherical' and 'cubic'"
-                    .format(selection.lower()))
-        elif isinstance(selection, np.array):
-            raise NotImplementedError(
-                "Numpy array indexing for subdomain extraction coming soon...")
-        else:
-            raise ValueError(
-                "Please specify a valid selection criterion (either 'spherical' or 'cubic')"
-            )
-
-        if not np.any(valid_comp_amr):
-            raise ValueError("Empty domain")
-
-        return valid_comp_amr, valid_comp_sinks, valid_comp_particles
-
-    def _extract_data(self, valid_comp_amr, valid_comp_sinks, valid_comp_particles):
-        # proceed to data extraction from parent dataset to subdomain
-        self.groups = {}
-        # extract amr dependent groups
-        for group in ["hydro", "amr", "grav"]:
-            self.groups[group] = Datagroup(parent=self, name=group)
-            for element in self._parent[group]:
-                self.groups[group][element] = self._parent[group][element][
-                    valid_comp_amr]
-        # extract sinks & particles
-        if self._sink_particles:
-            self.groups["sink"] = Datagroup(parent=self, name="sink")
-            for element in self._parent["sink"]:
-                self.groups["sink"][element] = self._parent["sink"][element][
-                    valid_comp_sinks]
-        if self._particles:
-            self.groups["part"] = Datagroup(parent=self, name="part")
-            for element in self._parent["part"]:
-                self.groups["part"][element] = self._parent["part"][element][
-                    valid_comp_particles]
diff --git a/src/osyris/io/sink.py b/src/osyris/io/sink.py
index 7d2b67f5..576fd7e1 100644
--- a/src/osyris/io/sink.py
+++ b/src/osyris/io/sink.py
@@ -50,7 +50,7 @@ def initialize(self, meta, units, select):
             else:
                 if all(x in u for x in ["[", "]"]):
                     # Legacy sink format quantities are not in code units
-                    unit_list.append(ureg(u.replace("[", "").replace("]", "")))
+                    unit_list.append(1.0 * ureg(u.replace("[", "").replace("]", "")))
                 else:
                     unit_list.append(eval(u.replace(' ', '*')))
 
diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py
index 2c699242..dc47ea67 100644
--- a/src/osyris/spatial/spatial.py
+++ b/src/osyris/spatial/spatial.py
@@ -4,172 +4,40 @@
 import numpy as np
 
 
-def change_origin(dataset, new_origin):
+def extract_sphere(dataset, radius, origin):
     """
-    Translate all positionnal coordinates to new origin
+    Extract a spherical subdomain around an origin point.
     """
-    for g in dataset.groups.keys():
-        for element in dataset[g]:
-            if "position" in element.lower():
-                dataset[g][element] -= new_origin
-    dataset.origin = new_origin
+    subdomain = dataset.__class__(nout=dataset.meta["nout"], path=dataset.meta["path"])
+    subdomain.meta = dataset.meta
+    subdomain._parent = dataset
 
+    for name, group in dataset.items():
+        pos = group.get("position", group.parent["amr"]["position"])
+        r = (pos - origin).norm
+        c = (r < radius).values
+        if np.any(c):
+            subdomain[name] = dataset[name][c]
 
-def rotation_matrix(vec, angle):
-    """
-    Returns 3D rotation matrix of angle 'angle' around rotation vector 'vec'.
-    """
-    if isinstance(vec, list):
-        vec = np.array(vec)
-    vec = vec / np.linalg.norm(vec)
-    R = np.cos(angle) * np.identity(3) + (np.sin(angle)) * np.cross(
-        vec,
-        np.identity(vec.shape[0]) * -1) + (1 - np.cos(angle)) * (np.outer(vec, vec))
-    return R
-
-
-def cartesian_to_spherical_rotation_matrix(colatitude, azimuth):
-    """
-    Used to convert arbitrary cartesian vectors into spherical vectors
-    """
-
-    R = np.array([[
-        np.sin(colatitude) * np.cos(azimuth),
-        np.sin(colatitude) * np.sin(azimuth),
-        np.cos(colatitude)
-    ],
-                  [
-                      np.cos(colatitude) * np.cos(azimuth),
-                      np.cos(colatitude) * np.sin(azimuth), -np.sin(colatitude)
-                  ], [-np.sin(azimuth), np.cos(azimuth), 0]],
-                 dtype=object)
-    return R
-
-
-def cartesian_to_cylindrical_rotation_matrix(azimuth):
-    """
-    Used to convert arbitrary cartesian vectors into cylindrical vectors
-    """
-
-    R = np.array([[np.cos(azimuth), np.sin(azimuth), 0],
-                  [-np.sin(azimuth), np.cos(azimuth), 0], [0, 0, 1]],
-                 dtype=object)
-    return R
-
-
-def get_ang_mom(subdomain, dr_L):
-    """
-    Compute angular momentum vector in sphere of radius dr_L
-    """
-    sphere = (subdomain["amr"]["position"].norm <= dr_L).values
-    pos = subdomain["amr"]["position"][sphere]
-    mv = subdomain["hydro"]["mass"][sphere] * subdomain["hydro"]["velocity"][sphere]
-    L = np.sum(pos.cross(mv))
-    return L
-
-
-def change_basis(subdomain, new_basis):
-    """
-    Rotates all vectors in dataset to align with vector in 'new_basis'
-    """
-    try:
-        old_basis = subdomain.basis
-    except AttributeError:
-        old_basis = [0, 0, 1]  # assume it's the ramses grid
-    if hasattr(new_basis, 'nvec'):  # if it's a vector
-        new_basis = [new_basis.x.values, new_basis.y.values, new_basis.z.values]
-    rot_angle = np.arccos(
-        np.dot(old_basis, new_basis) /
-        (np.linalg.norm(old_basis) * np.linalg.norm(new_basis)))
-    rot_vector = np.cross(new_basis, old_basis)
-    R = rotation_matrix(rot_vector, rot_angle)
-    for g in subdomain.groups.keys():
-        for element in subdomain[g]:
-            if hasattr(subdomain[g][element], 'nvec'):
-                # all of this will be simplified once matmul is integraded into Vector
-                vector = np.array([
-                    subdomain[g][element].x.values, subdomain[g][element].y.values,
-                    subdomain[g][element].z.values
-                ])
-                vector = R @ vector
-                subdomain[g][element] = subdomain[g][element].__class__(
-                    x=vector[0],
-                    y=vector[1],
-                    z=vector[2],
-                    unit=subdomain[g][element].unit)
-    subdomain.basis = new_basis
-
-
-def get_spherical_colatitude(position_vec):
-    """
-    Compute colatitude from cartesian position vector
-    """
-    X = position_vec.x.values
-    Y = position_vec.y.values
-    Z = position_vec.z.values
-
-    theta = np.arctan2(np.linalg.norm([X, Y], axis=0), Z)
-
-    return position_vec.x.__class__(values=theta, unit='rad', name="position_phi")
-
-
-def get_spherical_azimuth(position_vec):
-    """
-    Compute azimuth from cartesian position vector
-    """
-    X = position_vec.x.values
-    Y = position_vec.y.values
-    phi = np.arctan2(Y, X)
-
-    return position_vec.x.__class__(values=phi, unit='rad', name="position_phi")
-
-
-def get_cylindrical_radial_component(vec):
-    """
-    Get the radial component of an arbitrary cartesian vector 'vec' in cylindrical coords
-    """
-    if not hasattr(vec, 'nvec'):
-        raise ValueError("This is not a vector")
-    if vec.parent.name == "hydro":
-        position = vec.parent.parent["amr"]["position"]
-    else:
-        position = vec.parent["position"]
-
-    phi = position.phi
-    vec_r = np.cos(phi)*vec.x + np.sin(phi)*vec.y
-    vec_r.name = vec.name + "_cyl_r"
+    return subdomain
 
-    return vec_r
 
-def get_spherical_components(vec, comp):
+def extract_cube(dataset, dx, dy, dz, origin):
     """
-    Get the spherical components of an arbitrary cartesian vector 'vec'
+    Extract a cubic domain of size dx, dy & dz around an origin point
     """
-    if not hasattr(vec, 'nvec'):
-        raise ValueError("This is not a vector")
-    if vec.parent.name == "hydro":
-        position = vec.parent.parent["amr"]["position"]
-    else:
-        position = vec.parent["position"]
+    subdomain = dataset.__class__(nout=dataset.meta["nout"], path=dataset.meta["path"])
+    subdomain.meta = dataset.meta
+    subdomain._parent = dataset
 
-    theta = position.theta
-    phi = position.phi
+    for name, group in dataset.items():
+        pos = group.get("position", group.parent["amr"]["position"])
+        centered_pos = pos - origin
+        cx = (centered_pos.x <= dx * .5) & (centered_pos.x >= -dx * .5)
+        cy = (centered_pos.y <= dy * .5) & (centered_pos.y >= -dy * .5)
+        cz = (centered_pos.z <= dz * .5) & (centered_pos.z >= -dz * .5)
+        c = (cx & cy & cz).values
+        if np.any(c):
+            subdomain[name] = dataset[name][c]
 
-    if comp.lower() == "radius":
-        vec_r = np.cos(phi) * np.sin(theta) * vec.x + np.sin(phi) * np.sin(
-            theta) * vec.y + np.cos(theta) * vec.z
-        vec_r.name = vec.name + "_r"
-        return vec_r
-    elif comp.lower() == "colatitude":
-        vec_theta = np.cos(theta) * np.cos(phi) * vec.x + np.cos(theta) * np.sin(
-            phi) * vec.y - np.sin(theta) * vec.z
-        vec_theta.name = vec.name + "_theta"
-        return vec_theta
-    elif comp.lower() == "azimuth":
-        vec_phi = -np.sin(phi) * vec.x + np.cos(phi) * vec.y
-        vec_phi.name = vec.name + "_phi"
-        return vec_phi
-    else:
-        raise ValueError(
-            "Urecognized component keyword '{}'. Valid strings are 'radius', 'colatitude' and 'azimuth'"
-            .format(comp))
+    return subdomain
diff --git a/src/osyris/units/__init__.py b/src/osyris/units/__init__.py
new file mode 100644
index 00000000..5e65f58e
--- /dev/null
+++ b/src/osyris/units/__init__.py
@@ -0,0 +1,7 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
+
+# flake8: noqa
+
+from .units import units
+from .library import UnitsLibrary
diff --git a/src/osyris/units/library.py b/src/osyris/units/library.py
new file mode 100644
index 00000000..a6a7c052
--- /dev/null
+++ b/src/osyris/units/library.py
@@ -0,0 +1,45 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
+
+import re
+
+
+class UnitsLibrary:
+    """
+    A small helper class that holds the list of units, and behaves a bit like a
+    dict but also performs regex matching on key with wildcard character '*'.
+    """
+    def __init__(self, library, default_unit):
+        self._library = library
+        self._default_unit = default_unit
+
+    def __getitem__(self, key):
+        if key in self._library:
+            return self._library[key]
+        for name, unit in self._library.items():
+            if '*' in name:
+                regex = re.compile(name.replace('*', '.+'))
+                if re.match(regex, key):
+                    return unit
+        return self._default_unit
+
+    def __setitem__(self, key, value):
+        self._library[key] = value
+
+    def __str__(self):
+        return str(self._library)
+
+    def __repr__(self):
+        return str(self)
+
+    def keys(self):
+        return self._library.keys()
+
+    def values(self):
+        return self._library.values()
+
+    def items(self):
+        return self._library.items()
+
+    def update(self, *args, **kwargs):
+        self._library.update(*args, **kwargs)
diff --git a/src/osyris/units.py b/src/osyris/units/units.py
similarity index 96%
rename from src/osyris/units.py
rename to src/osyris/units/units.py
index 14ffd766..16ba4445 100644
--- a/src/osyris/units.py
+++ b/src/osyris/units/units.py
@@ -2,7 +2,7 @@
 # Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
 
 from pint import Quantity, UnitRegistry, Unit
-from . import config
+from .. import config
 
 
 class Units:
diff --git a/test/test_array.py b/test/test_array.py
index 7451759a..f92451b2 100644
--- a/test/test_array.py
+++ b/test/test_array.py
@@ -62,6 +62,15 @@ def test_constructor_masked_array():
     assert np.array_equal(array.values.mask, [False, False, False, True, True])
 
 
+def test_constructor_2d():
+    a = np.arange(12.).reshape(4, 3)
+    array = Array(values=a, unit='m')
+    assert array.unit == units('m')
+    assert len(array) == len(a)
+    assert array.shape == a.shape
+    assert np.array_equal(array.values, a)
+
+
 def test_addition():
     a = Array(values=[1., 2., 3., 4., 5.], unit='m')
     b = Array(values=[6., 7., 8., 9., 10.], unit='m')
@@ -636,22 +645,37 @@ def test_max():
 def test_reshape():
     a = Array(values=[1., 2., 3., 4., 5., 6.], unit='m')
     expected = Array(values=[[1., 2., 3.], [4., 5., 6.]], unit='m')
-    assert arraytrue(np.ravel(a.reshape(2, 3) == expected))
+    assert arrayequal(a.reshape(2, 3), expected)
 
 
 def test_slicing():
     a = Array(values=[11., 12., 13., 14., 15.], unit='m')
     assert a[2] == Array(values=[13.], unit='m')
-    assert arraytrue(a[:4] == Array(values=[11., 12., 13., 14.], unit='m'))
-    assert arraytrue(a[2:4] == Array(values=[13., 14.], unit='m'))
+    assert arrayequal(a[:4], Array(values=[11., 12., 13., 14.], unit='m'))
+    assert arrayequal(a[2:4], Array(values=[13., 14.], unit='m'))
 
 
-def test_slicing_vector():
+def test_slicing_2d():
     a = Array(values=np.arange(12.).reshape(4, 3), unit='m')
-    assert arraytrue(np.ravel(a[2:3] == Array(values=[[6., 7., 8.]], unit='m')))
+    assert arrayequal(a[0, :], Array(values=[0., 1., 2.], unit='m'))
+    assert arrayequal(a[1], Array(values=[3., 4., 5.], unit='m'))
+    assert arrayequal(a[2:3], Array(values=[[6., 7., 8.]], unit='m'))
     assert a[2:3].shape == (1, 3)
-    assert arraytrue(
-        np.ravel(a[:2] == Array(values=[[0., 1., 2.], [3., 4., 5.]], unit='m')))
+    assert arrayequal(a[:2], Array(values=[[0., 1., 2.], [3., 4., 5.]], unit='m'))
+
+
+def test_slicing_with_array():
+    a = Array(values=[1., 2., 3., 4., 5.], unit='m')
+    select = a > (3. * units('m'))
+    assert arrayequal(a[select], Array(values=[4., 5.], unit='m'))
+    b = Array(values=[0, 2, 4])
+    assert arrayequal(a[b], Array(values=[1., 3., 5.], unit='m'))
+
+
+def test_slicing_with_array_bad_dtype():
+    a = Array(values=[1., 2., 3., 4., 5.], unit='m')
+    with pytest.raises(TypeError):
+        _ = a[Array(values=[0., 2., 4.])]
 
 
 def test_copy():
diff --git a/test/test_vector.py b/test/test_vector.py
index c7545195..1e8dd5ab 100644
--- a/test/test_vector.py
+++ b/test/test_vector.py
@@ -713,6 +713,29 @@ def test_slicing():
     assert vectorequal(v[2:4], Vector(x=x[2:4], y=y[2:4], z=z[2:4]))
 
 
+def test_slicing_with_array():
+    x = Array(values=[1., 2., 3., 4., 5.], unit='m')
+    y = Array(values=[6., 7., 8., 9., 10.], unit='m')
+    z = Array(values=[11., 12., 13., 14., 15.], unit='m')
+    v = Vector(x, y, z)
+    select = v.norm > (15. * units('m'))
+    exp_x = Array(values=[3., 4., 5.], unit='m')
+    exp_y = Array(values=[8., 9., 10.], unit='m')
+    exp_z = Array(values=[13., 14., 15.], unit='m')
+    expected = Vector(exp_x, exp_y, exp_z)
+    assert vectorequal(v[select], expected)
+
+
+def test_slicing_with_vector_raises():
+    x = Array(values=[1., 2., 3., 4., 5.], unit='m')
+    y = Array(values=[6., 7., 8., 9., 10.], unit='m')
+    z = Array(values=[11., 12., 13., 14., 15.], unit='m')
+    v = Vector(x, y, z)
+    select = v > (16. * units('m'))
+    with pytest.raises(ValueError):
+        _ = v[select]
+
+
 def test_copy():
     x = Array(values=[1., 2., 3., 4., 5.], unit='m')
     y = Array(values=[6., 7., 8., 9., 10.], unit='m')

From e2018f05d11559728b2a43ecedc441c613950f8e Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 13 Jun 2022 16:58:13 +0200
Subject: [PATCH 14/28] add new line

---
 src/osyris/core/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/osyris/core/__init__.py b/src/osyris/core/__init__.py
index e5ab547a..bb26aeca 100644
--- a/src/osyris/core/__init__.py
+++ b/src/osyris/core/__init__.py
@@ -7,4 +7,4 @@
 from .vector import Vector
 from .datagroup import Datagroup
 from .dataset import Dataset
-from .plot import Plot
\ No newline at end of file
+from .plot import Plot

From 3af99907f7f2faf94086f76dfa57735c2ce91003 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 13 Jun 2022 16:58:26 +0200
Subject: [PATCH 15/28] revert name assignement change

---
 src/osyris/core/datagroup.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/osyris/core/datagroup.py b/src/osyris/core/datagroup.py
index 5a0f882f..421ca597 100644
--- a/src/osyris/core/datagroup.py
+++ b/src/osyris/core/datagroup.py
@@ -5,10 +5,10 @@
 
 
 class Datagroup:
-    def __init__(self, data=None, parent=None, name=""):
+    def __init__(self, data=None, parent=None):
         self._container = {}
         self._parent = parent
-        self._name = name
+        self._name = ""
         self.shape = None
         if data is not None:
             for key, array in data.items():

From 015d0064afe88c9baad10f594fd0ead056cbcfe3 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 13 Jun 2022 16:58:43 +0200
Subject: [PATCH 16/28] revert spherical component computations for now

---
 src/osyris/core/vector.py | 33 ---------------------------------
 1 file changed, 33 deletions(-)

diff --git a/src/osyris/core/vector.py b/src/osyris/core/vector.py
index f0f46acd..ad10ddf7 100644
--- a/src/osyris/core/vector.py
+++ b/src/osyris/core/vector.py
@@ -5,7 +5,6 @@
 from .base import Base
 from .array import Array
 from .tools import value_to_string
-from .. import spatial
 
 
 def _binary_op(op, lhs, rhs):
@@ -254,35 +253,3 @@ def cross(self, other):
         z = self.x * other.y
         z -= self.y * other.x
         return self.__class__(x, y, z)
-
-    @property
-    def r(self):
-        if self.name == "position":
-            v = self.norm
-            v.name = "position_r"
-            return v
-        else:
-            return spatial.get_spherical_components(self, comp='radius')
-
-    @property
-    def theta(self):
-        if self.name == "position":
-            return spatial.get_spherical_colatitude(self)
-        else:
-            return spatial.get_spherical_components(self, comp='colatitude')
-
-    @property
-    def phi(self):
-        if self.name == "position":
-            return spatial.get_spherical_azimuth(self)
-        else:
-            return spatial.get_spherical_components(self, comp='azimuth')
-
-    @property
-    def cyl_r(self):
-        if self.name == "position":
-            v = Array(values=np.sqrt(self.x.values**2 + self.y.values**2), unit=self.unit)
-            v.name = "position_cyl_r"
-            return v
-        else:
-            return spatial.get_cylindrical_radial_component(self)
\ No newline at end of file

From 8faa3ffe697a51027745ac4c368abd48343a27ae Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 13 Jun 2022 17:09:52 +0200
Subject: [PATCH 17/28] newline

---
 src/osyris/spatial/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py
index 5b6a590a..307f08a5 100644
--- a/src/osyris/spatial/__init__.py
+++ b/src/osyris/spatial/__init__.py
@@ -3,4 +3,4 @@
 
 # flake8: noqa
 
-from .spatial import *
\ No newline at end of file
+from .spatial import *

From 81ef3e5bd00abc125b395fb04e2a5c20409d9a0a Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <60321972+Adnan-Ali-Ahmad@users.noreply.github.com>
Date: Tue, 14 Jun 2022 16:46:59 +0200
Subject: [PATCH 18/28] Update src/osyris/spatial/__init__.py

Co-authored-by: Neil Vaytet <39047984+nvaytet@users.noreply.github.com>
---
 src/osyris/spatial/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py
index 307f08a5..d0c17376 100644
--- a/src/osyris/spatial/__init__.py
+++ b/src/osyris/spatial/__init__.py
@@ -1,5 +1,5 @@
 # SPDX-License-Identifier: BSD-3-Clause
-# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris)
+# Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
 
 # flake8: noqa
 

From 4be4096662013768f3f61e9146756acc592a2178 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <60321972+Adnan-Ali-Ahmad@users.noreply.github.com>
Date: Tue, 14 Jun 2022 16:47:20 +0200
Subject: [PATCH 19/28] Update src/osyris/spatial/spatial.py

Co-authored-by: Neil Vaytet <39047984+nvaytet@users.noreply.github.com>
---
 src/osyris/spatial/spatial.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py
index dc47ea67..e16d7a22 100644
--- a/src/osyris/spatial/spatial.py
+++ b/src/osyris/spatial/spatial.py
@@ -1,5 +1,5 @@
 # SPDX-License-Identifier: BSD-3-Clause
-# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris)
+# Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
 
 import numpy as np
 

From a88c647253142d7c97d5e5b212bfe721ec8cb601 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <60321972+Adnan-Ali-Ahmad@users.noreply.github.com>
Date: Tue, 14 Jun 2022 16:49:40 +0200
Subject: [PATCH 20/28] Update src/osyris/spatial/spatial.py

Co-authored-by: Neil Vaytet <39047984+nvaytet@users.noreply.github.com>
---
 src/osyris/spatial/spatial.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py
index e16d7a22..6c86e96a 100644
--- a/src/osyris/spatial/spatial.py
+++ b/src/osyris/spatial/spatial.py
@@ -22,7 +22,7 @@ def extract_sphere(dataset, radius, origin):
     return subdomain
 
 
-def extract_cube(dataset, dx, dy, dz, origin):
+def extract_box(dataset, dx, dy, dz, origin):
     """
     Extract a cubic domain of size dx, dy & dz around an origin point
     """

From 2d6b80554ebcdc18e52871435942f3989a898bc2 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Tue, 14 Jun 2022 17:01:25 +0200
Subject: [PATCH 21/28] renaming spatial.py to subdomain.py

---
 src/osyris/spatial/__init__.py                  | 2 +-
 src/osyris/spatial/{spatial.py => subdomain.py} | 0
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename src/osyris/spatial/{spatial.py => subdomain.py} (100%)

diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py
index d0c17376..4f59a55f 100644
--- a/src/osyris/spatial/__init__.py
+++ b/src/osyris/spatial/__init__.py
@@ -3,4 +3,4 @@
 
 # flake8: noqa
 
-from .spatial import *
+from .subdomain import extract_sphere, extract_box
diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/subdomain.py
similarity index 100%
rename from src/osyris/spatial/spatial.py
rename to src/osyris/spatial/subdomain.py

From b5c3d6eaf09f44a606fb3fc013731a0495d493e4 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Tue, 14 Jun 2022 17:04:33 +0200
Subject: [PATCH 22/28] reverted name setter in .to method

---
 src/osyris/core/vector.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/osyris/core/vector.py b/src/osyris/core/vector.py
index ad10ddf7..a37824c9 100644
--- a/src/osyris/core/vector.py
+++ b/src/osyris/core/vector.py
@@ -209,9 +209,7 @@ def __invert__(self):
         return np.logical_not(self)
 
     def to(self, unit):
-        return self.__class__(**{c: xyz.to(unit)
-                                 for c, xyz in self._xyz.items()},
-                              name=self.name)
+        return self.__class__(**{c: xyz.to(unit) for c, xyz in self._xyz.items()})
 
     def _wrap_numpy(self, func, *args, **kwargs):
         if isinstance(args[0], (tuple, list)):

From 668789af5a412884707b1432ba122a2ea36b21a5 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Wed, 15 Jun 2022 15:39:32 +0200
Subject: [PATCH 23/28] removing parent and meta attribution

---
 src/osyris/__init__.py          | 1 +
 src/osyris/spatial/subdomain.py | 4 ----
 2 files changed, 1 insertion(+), 4 deletions(-)

diff --git a/src/osyris/__init__.py b/src/osyris/__init__.py
index 67353306..25114bb3 100644
--- a/src/osyris/__init__.py
+++ b/src/osyris/__init__.py
@@ -7,3 +7,4 @@
 from .units import units
 from .core import Array, Datagroup, Dataset, Plot, Vector
 from .plot import histogram1d, histogram2d, plane, scatter, map, plot
+from .spatial import extract_box, extract_sphere
diff --git a/src/osyris/spatial/subdomain.py b/src/osyris/spatial/subdomain.py
index 6c86e96a..c2c4b47c 100644
--- a/src/osyris/spatial/subdomain.py
+++ b/src/osyris/spatial/subdomain.py
@@ -9,8 +9,6 @@ def extract_sphere(dataset, radius, origin):
     Extract a spherical subdomain around an origin point.
     """
     subdomain = dataset.__class__(nout=dataset.meta["nout"], path=dataset.meta["path"])
-    subdomain.meta = dataset.meta
-    subdomain._parent = dataset
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])
@@ -27,8 +25,6 @@ def extract_box(dataset, dx, dy, dz, origin):
     Extract a cubic domain of size dx, dy & dz around an origin point
     """
     subdomain = dataset.__class__(nout=dataset.meta["nout"], path=dataset.meta["path"])
-    subdomain.meta = dataset.meta
-    subdomain._parent = dataset
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])

From cd909e760c462eb3cf23c8eba38416bcf4f5aa20 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Fri, 24 Jun 2022 14:09:51 +0200
Subject: [PATCH 24/28] removing loader from subdomain

---
 src/osyris/spatial/subdomain.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/osyris/spatial/subdomain.py b/src/osyris/spatial/subdomain.py
index c2c4b47c..9cf466d4 100644
--- a/src/osyris/spatial/subdomain.py
+++ b/src/osyris/spatial/subdomain.py
@@ -8,7 +8,8 @@ def extract_sphere(dataset, radius, origin):
     """
     Extract a spherical subdomain around an origin point.
     """
-    subdomain = dataset.__class__(nout=dataset.meta["nout"], path=dataset.meta["path"])
+    subdomain = dataset.__class__()
+    subdomain.meta = dataset.meta.copy()
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])
@@ -24,7 +25,8 @@ def extract_box(dataset, dx, dy, dz, origin):
     """
     Extract a cubic domain of size dx, dy & dz around an origin point
     """
-    subdomain = dataset.__class__(nout=dataset.meta["nout"], path=dataset.meta["path"])
+    subdomain = dataset.__class__()
+    subdomain.meta = dataset.meta.copy()
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])

From 6fd890372ac031200c1193d25b4bec7be522f6f0 Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 27 Jun 2022 08:35:49 +0200
Subject: [PATCH 25/28] implemented shape check

---
 src/osyris/spatial/subdomain.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/osyris/spatial/subdomain.py b/src/osyris/spatial/subdomain.py
index 9cf466d4..af926b07 100644
--- a/src/osyris/spatial/subdomain.py
+++ b/src/osyris/spatial/subdomain.py
@@ -13,6 +13,8 @@ def extract_sphere(dataset, radius, origin):
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])
+        if pos.shape[0] != group.shape:
+            continue
         r = (pos - origin).norm
         c = (r < radius).values
         if np.any(c):
@@ -30,6 +32,8 @@ def extract_box(dataset, dx, dy, dz, origin):
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])
+        if pos.shape[0] != group.shape:
+            continue
         centered_pos = pos - origin
         cx = (centered_pos.x <= dx * .5) & (centered_pos.x >= -dx * .5)
         cy = (centered_pos.y <= dy * .5) & (centered_pos.y >= -dy * .5)

From 98b87021d69df6dabac223d0d3dfe8fb64ad919d Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Mon, 27 Jun 2022 08:39:30 +0200
Subject: [PATCH 26/28] import dataset

---
 src/osyris/spatial/subdomain.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/osyris/spatial/subdomain.py b/src/osyris/spatial/subdomain.py
index af926b07..a45cbf19 100644
--- a/src/osyris/spatial/subdomain.py
+++ b/src/osyris/spatial/subdomain.py
@@ -2,13 +2,14 @@
 # Copyright (c) 2022 Osyris contributors (https://github.com/osyris-project/osyris)
 
 import numpy as np
+from .. import Dataset
 
 
 def extract_sphere(dataset, radius, origin):
     """
     Extract a spherical subdomain around an origin point.
     """
-    subdomain = dataset.__class__()
+    subdomain = Dataset()
     subdomain.meta = dataset.meta.copy()
 
     for name, group in dataset.items():
@@ -27,7 +28,7 @@ def extract_box(dataset, dx, dy, dz, origin):
     """
     Extract a cubic domain of size dx, dy & dz around an origin point
     """
-    subdomain = dataset.__class__()
+    subdomain = Dataset()
     subdomain.meta = dataset.meta.copy()
 
     for name, group in dataset.items():

From 16e522a47a91d9d5429e19f24fe7ca90b96cfdeb Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Thu, 30 Jun 2022 13:45:01 +0200
Subject: [PATCH 27/28] warning message added

---
 src/osyris/spatial/subdomain.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/osyris/spatial/subdomain.py b/src/osyris/spatial/subdomain.py
index a45cbf19..ad4e47a9 100644
--- a/src/osyris/spatial/subdomain.py
+++ b/src/osyris/spatial/subdomain.py
@@ -3,6 +3,7 @@
 
 import numpy as np
 from .. import Dataset
+import warnings
 
 
 def extract_sphere(dataset, radius, origin):
@@ -14,7 +15,10 @@ def extract_sphere(dataset, radius, origin):
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])
-        if pos.shape[0] != group.shape:
+        if pos.shape != group.shape:
+            warnings.warn(
+                "Ignoring datagroup '{}', which has no position " +
+                "vector and has different shape than 'amr' group.".format(group))
             continue
         r = (pos - origin).norm
         c = (r < radius).values
@@ -33,7 +37,10 @@ def extract_box(dataset, dx, dy, dz, origin):
 
     for name, group in dataset.items():
         pos = group.get("position", group.parent["amr"]["position"])
-        if pos.shape[0] != group.shape:
+        if pos.shape != group.shape:
+            warnings.warn(
+                "Ignoring datagroup '{}', which has no position " +
+                "vector and has different shape than 'amr' group.".format(group))
             continue
         centered_pos = pos - origin
         cx = (centered_pos.x <= dx * .5) & (centered_pos.x >= -dx * .5)

From 398962c4f571c0545bf61f4b013d5074177daadd Mon Sep 17 00:00:00 2001
From: Adnan-Ali-Ahmad <adnan.ali.ahmad1998@gmail.com>
Date: Thu, 30 Jun 2022 13:56:38 +0200
Subject: [PATCH 28/28] flake8

---
 src/osyris/spatial/subdomain.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/osyris/spatial/subdomain.py b/src/osyris/spatial/subdomain.py
index ad4e47a9..3ea90035 100644
--- a/src/osyris/spatial/subdomain.py
+++ b/src/osyris/spatial/subdomain.py
@@ -17,8 +17,8 @@ def extract_sphere(dataset, radius, origin):
         pos = group.get("position", group.parent["amr"]["position"])
         if pos.shape != group.shape:
             warnings.warn(
-                "Ignoring datagroup '{}', which has no position " +
-                "vector and has different shape than 'amr' group.".format(group))
+                "Ignoring datagroup '{}', which has no position ".format(group) +
+                "vector and has different shape than 'amr' group.")
             continue
         r = (pos - origin).norm
         c = (r < radius).values
@@ -39,8 +39,8 @@ def extract_box(dataset, dx, dy, dz, origin):
         pos = group.get("position", group.parent["amr"]["position"])
         if pos.shape != group.shape:
             warnings.warn(
-                "Ignoring datagroup '{}', which has no position " +
-                "vector and has different shape than 'amr' group.".format(group))
+                "Ignoring datagroup '{}', which has no position ".format(group) +
+                "vector and has different shape than 'amr' group.")
             continue
         centered_pos = pos - origin
         cx = (centered_pos.x <= dx * .5) & (centered_pos.x >= -dx * .5)