Merge pull request #78 from kyleaoman/fix_tng

Fix a bug in nH_g calculation for TNG.
kyleaoman · Oct 12, 2024 · 5fb378d · 5fb378d
2 parents b12db96 + 77e7265
commit 5fb378d
Show file tree

Hide file tree

Showing 10 changed files with 1,853 additions and 347 deletions.
diff --git a/.gitignore b/.gitignore
@@ -18,4 +18,5 @@ examples/*.tar
 examples/sql_cookies.txt
 examples/eagle_db.password
 examples/tng_api.key
+examples/m11e_res7100
 .coverage
diff --git a/codemeta.json b/codemeta.json
@@ -19,6 +19,6 @@
     "codeRepository": [
         "https://github.com/kyleaoman/martini",
     ],
-    "version": "2.1.5",
+    "version": "2.1.6",
     "license": "https://spdx.org/licenses/GPL-3.0-only.html",
 }
diff --git a/examples/martini_TNG.ipynb b/examples/martini_TNG.ipynb
diff --git a/examples/martini_eagle.ipynb b/examples/martini_eagle.ipynb
diff --git a/examples/martini_fire.ipynb b/examples/martini_fire.ipynb
diff --git a/examples/martini_simba.ipynb b/examples/martini_simba.ipynb
diff --git a/examples/martini_source_injection.ipynb b/examples/martini_source_injection.ipynb
diff --git a/martini/__version__.py b/martini/__version__.py
@@ -1 +1 @@
-__version__ = "2.1.5"
+__version__ = "2.1.6"
diff --git a/martini/sources/tng_source.py b/martini/sources/tng_source.py
@@ -348,17 +348,32 @@ def __init__(
         xe_g = data_g["ElectronAbundance"]
         rho_g = data_g["Density"] * 1e10 / h * U.Msun * np.power(a / h * U.kpc, -3)
         u_g = data_g["InternalEnergy"]  # unit conversion handled in T_g
-        mu_g = 4 * C.m_p.to(U.g).value / (1 + 3 * X_H_g + 4 * X_H_g * xe_g)
+        mu_g = 4 / (1 + 3 * X_H_g + 4 * X_H_g * xe_g)
         gamma = 5.0 / 3.0  # see http://www.tng-project.org/data/docs/faq/#gen4
-        T_g = (gamma - 1) * u_g / C.k_B.to(U.erg / U.K).value * 1e10 * mu_g * U.K
+        T_g = (
+            (gamma - 1)
+            * u_g
+            / C.k_B.to_value(U.erg / U.K)
+            * 1e10
+            * mu_g
+            * C.m_p.to_value(U.g)
+            * U.K
+        )
         m_g = data_g["Masses"] * 1e10 / h * U.Msun
         # cast to float64 to avoid underflow error
         nH_g = U.Quantity(rho_g * X_H_g / mu_g, dtype=np.float64) / C.m_p
         # In TNG_corrections I set f_neutral = 1 for particles with density
         # > .1cm^-3. Might be possible to do a bit better here, but HI & H2
         # tables for TNG will be available soon anyway.
         fatomic_g = atomic_frac(
-            z, nH_g, T_g, rho_g, X_H_g, onlyA1=True, TNG_corrections=True
+            z,
+            nH_g,
+            T_g,
+            rho_g,
+            X_H_g,
+            mu=mu_g,
+            onlyA1=True,
+            TNG_corrections=True,
         )
         mHI_g = m_g * X_H_g * fatomic_g
         try:

diff --git a/pyproject.toml b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "astromartini"
-version = "2.1.5"
+version = "2.1.6"
 authors = [
     { name="Kyle Oman", email="[email protected]" },
 ]