clusters/help: Add radius_stripping()

This function calculates the stripping radius of the sub-cluster when it
falling through the main cluster.

1 files changed, 73 insertions, 6 deletions

diff --git a/fg21sim/extragalactic/clusters/helper.py b/fg21sim/extragalactic/clusters/helper.py
index b58dc87..ecbb866 100644
--- a/fg21sim/extragalactic/clusters/helper.py
+++ b/fg21sim/extragalactic/clusters/helper.py
@@ -51,6 +51,7 @@ import logging
 
 import numpy as np
 from scipy import integrate
+from scipy import optimize
 
 from ...share import CONFIGS, COSMO
 from ...utils.units import (Units as AU,
@@ -63,6 +64,46 @@ from ...utils.transform import circle2ellipse
 logger = logging.getLogger(__name__)
 
 
+def beta_model(rho0, rc, beta):
+    """
+    Return a function that calculates the value (gas density) at a radius
+    according to the β-model with the given parameters.
+    """
+    def func(r):
+        x = r / rc
+        return rho0 * (1 + x**2) ** (-1.5*beta)
+
+    return func
+
+
+def calc_gas_density_profile(mass, z):
+    """
+    Calculate the parameters of the β-model that is used to describe the
+    gas density profile.
+
+    NOTE
+    ----
+    The core radius is assumed to be: ``rc = 0.1 r_vir``, and the beta
+    parameters is assumed to be ``β = 0.8``.
+
+    Reference: [cassano2005],Sec.(4.1)
+
+    Returns
+    -------
+    fbeta : function
+        A function of the β-model.
+        Unit: [Msun/kpc^3]
+    """
+    r_vir = radius_virial(mass, z)  # [kpc]
+    rc = 0.1 * r_vir
+    beta = 0.8
+    fint = beta_model(1, rc, beta)
+    v = integrate.quad(lambda r: fint(r) * r**2,
+                       a=0, b=r_vir)[0]  # [kpc^3]
+    rho0 = mass * COSMO.baryon_fraction / (4*np.pi * v)  # [Msun/kpc^3]
+    return beta_model(rho0, rc, beta)
+
+
 def radius_virial(mass, z=0.0):
     """
     Calculate the virial radius of a cluster at a given redshift.
@@ -85,8 +126,7 @@ def radius_virial(mass, z=0.0):
     Dc = COSMO.overdensity_virial(z)
     rho = COSMO.rho_crit(z)  # [g/cm^3]
     R_vir = (3*mass*AUC.Msun2g / (4*np.pi * Dc * rho))**(1/3)  # [cm]
-    R_vir *= AUC.cm2kpc  # [kpc]
-    return R_vir
+    return R_vir * AUC.cm2kpc  # [kpc]
 
 
 def radius_halo(mass, z=0.0, configs=CONFIGS):
@@ -124,6 +164,35 @@ def radius_halo(mass, z=0.0, configs=CONFIGS):
     return R_halo
 
 
+def radius_stripping(M_main, M_sub, z, configs=CONFIGS):
+    """
+    Calculate the stripping radius of the in-falling sub-cluster, which
+    is determined by the equipartition between the static and ram pressure.
+
+    Reference: [cassano2005],Eqs.(11,12,13,14)
+
+    Returns
+    -------
+    rs : float
+        The stripping radius of the sub-cluster.
+        Unit: [kpc]
+    """
+    r_vir = radius_virial(M_sub, z)  # [kpc]
+    rho_main = density_number_thermal(M_main, z) * AC.mu*AC.u  # [g/cm^3]
+    f_rho_sub = calc_gas_density_profile(M_sub, z)  # [Msun/kpc^3]
+    vi = velocity_impact(M_main, M_sub, z)  # [km/s]
+    kT_sub = kT_cluster(M_sub, z)  # [keV]
+    rhs = rho_main * vi**2 * AC.mu*AC.u / kT_sub  # [g/cm^3][g*km^2/s^2/keV]
+    rhs *= 1e3 * AUC.J2erg / AUC.keV2erg  # [g/cm^3]
+    rhs *= AUC.g2Msun / AUC.cm2kpc**3  # [Msun/kpc^3]
+    try:
+        rs = optimize.brentq(lambda r: f_rho_sub(r) - rhs,
+                             a=0.1*r_vir, b=r_vir, xtol=1e-1)
+    except ValueError:
+        rs = 0.1 * r_vir
+    return rs  # [kpc]
+
+
 def kT_virial(mass, z=0.0, radius=None):
     """
     Calculate the virial temperature of a cluster.
@@ -308,8 +377,7 @@ def velocity_virial(mass, z=0.0):
     """
     R_vir = radius_virial(mass, z) * AUC.kpc2cm  # [cm]
     vv = np.sqrt(AC.G * mass*AUC.Msun2g / R_vir)  # [cm/s]
-    vv /= AUC.km2cm  # [km/s]
-    return vv
+    return vv / AUC.km2cm  # [km/s]
 
 
 def velocity_impact(M_main, M_sub, z=0.0):
@@ -339,8 +407,7 @@ def velocity_impact(M_main, M_sub, z=0.0):
     R_vir = radius_virial(M_main, z) * AUC.kpc2cm  # [cm]
     vi = np.sqrt(2*AC.G * (1-1/eta_v) *
                  (M_main+M_sub)*AUC.Msun2g / R_vir)  # [cm/s]
-    vi /= AUC.km2cm  # [km/s]
-    return vi
+    return vi / AUC.km2cm  # [km/s]
 
 
 def time_crossing(M_main, M_sub, z=0.0):