From 1441e6441887b5f7e4b16801ced49aa65b245511 Mon Sep 17 00:00:00 2001
From: hfoote <haydenfoote@arizona.edu>
Date: Mon, 18 Nov 2024 11:43:00 -0700
Subject: [PATCH 1/2] Adding utils.units

---
 EXPtools/utils/units.py | 68 +++++++++++++++++++++++++++++++++++++++++
 tests/test_units.py     | 16 ++++++++++
 2 files changed, 84 insertions(+)
 create mode 100644 EXPtools/utils/units.py
 create mode 100644 tests/test_units.py

diff --git a/EXPtools/utils/units.py b/EXPtools/utils/units.py
new file mode 100644
index 0000000..17a72df
--- /dev/null
+++ b/EXPtools/utils/units.py
@@ -0,0 +1,68 @@
+"""
+Routines for defining and handling unit systems.
+"""
+
+import astropy.units as u
+import astropy.constants as const
+import numpy as np
+
+def define_exp_units(M_phys, R_phys, G, verbose=False):
+    """
+    Defines EXP (G=1) units in astropy for a given halo,
+    when given mass and length in physical units. 
+    Base units are called expM, expL, expV, expT. 
+
+    Parameters
+    ----------
+    M_phys : astropy quantity
+        Desired EXP mass unit in some physical unit of mass, 
+        i.e. a virial mass in solar masses or the mass unit of an external simulation
+        
+    R_phys : astropy quantity
+        Desired EXP length unit in some physical unit of length, 
+        i.e. a virial radius in kpc or the length unit of an external simulation
+        
+    G : astropy quantity
+        G in the physical unit system of the external simulation
+        
+    verbose : bool, optional
+        If True, print out the EXP unit system's conversion to physical units
+
+    Returns
+    -------
+    expL : astropy unit
+        EXP length unit
+        
+    expM : astropy unit
+        EXP mass unit
+        
+    expV : astropy unit
+        EXP velocity unit
+        
+    expT : astropy unit
+        EXP time unit
+    
+    Example
+    -------
+    To define the conversion from e.g. Gadget default units, use
+    expL, expM, expV, expT = define_exp_units(M_phys=1.*u.Msun, R_phys=1.*u.kpc, G=const.G.to(u.kpc/(1e10*u.Msun)*u.km**2/u.s**2))
+    """
+    
+    if verbose:
+        print("Using physical G = ", G)
+    
+    # define EXP velocity and time 
+    expL = u.def_unit('expL', R_phys)
+    expM = u.def_unit('expM', M_phys)
+    expV = u.def_unit('expV', np.sqrt(G*expM/expL))
+    expT = u.def_unit('expT', expL/expV)
+    
+    # print and return the unit system
+    if verbose:
+        print('Conversions between EXP units and physical units:')
+        print('1 expL = ', 1.*expL.to(u.kpc), 'kpc')
+        print('1 expM = ', 1.*expM.to(u.Msun), 'Msun')
+        print('1 expV = ', 1.*expV.to(u.km/u.s), 'km/s')
+        print('1 expT = ', 1.*expT.to(u.Gyr), 'Gyr')
+
+    return expL, expM, expV, expT
\ No newline at end of file
diff --git a/tests/test_units.py b/tests/test_units.py
new file mode 100644
index 0000000..f5ccc82
--- /dev/null
+++ b/tests/test_units.py
@@ -0,0 +1,16 @@
+"""
+Tests for unit conversions
+"""
+
+import numpy as np
+import astropy.units as u
+import astropy.constants as const
+from EXPtools.utils.exp_units import define_exp_units
+
+def test_define_exp_units():
+
+	tol = 1e-10
+	expL, expM, expV, expT = define_exp_units(M_phys=1.*u.Msun, R_phys=1.*u.kpc, G=const.G.to(u.kpc/(1e10*u.Msun)*u.km**2/u.s**2))
+	assert (np.abs(expT.value - (u.kpc/(u.km/u.s)).to(u.Gyr)) < tol), "Unit definitions in define_exp_units failed!"
+
+	return None

From 6c302b4fc04b10825a5964608b974e2162fa1dad Mon Sep 17 00:00:00 2001
From: Hayden Foote <haydenfoote@email.arizona.edu>
Date: Mon, 18 Nov 2024 12:08:18 -0700
Subject: [PATCH 2/2] fixed test for unit handling

---
 EXPtools/utils/__init__.py | 1 +
 tests/test_units.py        | 6 +++---
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/EXPtools/utils/__init__.py b/EXPtools/utils/__init__.py
index 546006b..a54c207 100644
--- a/EXPtools/utils/__init__.py
+++ b/EXPtools/utils/__init__.py
@@ -1,3 +1,4 @@
 from . import halo
 from . import coefficients
 from . import indexing
+from . import units
diff --git a/tests/test_units.py b/tests/test_units.py
index f5ccc82..966234d 100644
--- a/tests/test_units.py
+++ b/tests/test_units.py
@@ -5,12 +5,12 @@
 import numpy as np
 import astropy.units as u
 import astropy.constants as const
-from EXPtools.utils.exp_units import define_exp_units
+from EXPtools.utils import units
 
 def test_define_exp_units():
 
 	tol = 1e-10
-	expL, expM, expV, expT = define_exp_units(M_phys=1.*u.Msun, R_phys=1.*u.kpc, G=const.G.to(u.kpc/(1e10*u.Msun)*u.km**2/u.s**2))
-	assert (np.abs(expT.value - (u.kpc/(u.km/u.s)).to(u.Gyr)) < tol), "Unit definitions in define_exp_units failed!"
+	expL, expM, expV, expT = units.define_exp_units(M_phys=1.*u.Msun, R_phys=1.*u.kpc, G=const.G.to(u.kpc/(1e10*u.Msun)*u.km**2/u.s**2))
+	assert (np.abs((1.*expV).to(u.km/u.s).value - 0.002073865296984421) < tol), "Unit definitions in define_exp_units failed!"
 
 	return None