1 files changed, 71 insertions, 0 deletions
diff --git a/fg21sim/utils/virial.py b/fg21sim/utils/virial.py
new file mode 100644
index 0000000..2e81ab4
--- /dev/null
+++ b/fg21sim/utils/virial.py
@@ -0,0 +1,71 @@
+# Copyright (c) 2019 Weitian LI <wt@liwt.net>
+# MIT License
+
+"""
+Virial radius/mass utilities.
+
+.. [duffy2008]
+   Duffy et al. 2008, MNRAS, 390, L64
+   http://adsabs.harvard.edu/abs/2008MNRAS.390L..64D
+
+.. [ettori2009]
+   Ettori & Balestra 2009, A&A, 496, 343
+   http://adsabs.harvard.edu/abs/2009A%26A...496..343
+
+.. [lokas2001]
+   Łokas & Mamon 2001, MNRAS, 321, 155
+   http://adsabs.harvard.edu/abs/2001MNRAS.321..155L
+"""
+
+import logging
+
+import numpy as np
+from ..share import COSMO
+
+
+logger = logging.getLogger(__name__)
+
+
+def concentration(mass, z=0):
+    """
+    Calculate the NFW concentration parameter (c = r_vir / r_s) from the
+    mass-concentration relation: c = A [M/M_pivot]^B (1+z)^C
+
+    Reference: [duffy2008],Tab.1
+    """
+    A, B, C = 5.71, -0.084, -0.47  # for M200
+    M_pivot = 2e12 / COSMO.h  # [Msun]
+    return A * (mass/M_pivot) ** B + (1+z) ** C
+
+
+def mass_fraction_nfw(s, c):
+    """
+    Calculate the mass fraction (of virial mass) within the scaled radius
+    ``s = r / r_vir``.
+
+    Reference: [lokas2001]
+
+    Parameters
+    ----------
+    s : float
+        The radius scaled by the virial radius, i.e., s = r / r_vir.
+    c : float
+        The concentration parameter, i.e., c = r_vir / r_s.
+    """
+    nom = np.log(1 + c*s) - c*s / (1 + c*s)
+    denom = np.log(1 + c) - c / (1 + c)
+    return nom / denom
+
+
+def M200(M500, z=0):
+    """
+    Estimate the M200 from the M500.
+    """
+    # r500 ~= 0.65 r200; Ref.[ettori2009],Sec.2
+    s = 0.65
+    # Begin with 'M200 = 1.5*M500'
+    c = concentration(M500*1.5, z)
+    mass = M500 / mass_fraction_nfw(s, c)
+    # Again with the updated mass
+    c = concentration(mass, z)
+    return M500 / mass_fraction_nfw(s, c)