diff --git a/frb/halos/models.py b/frb/halos/models.py index 852a5a15..72165c33 100644 --- a/frb/halos/models.py +++ b/frb/halos/models.py @@ -1,7 +1,5 @@ """ Module for DM Halo calculations """ -from __future__ import print_function, absolute_import, division, unicode_literals - import numpy as np import pdb import warnings @@ -544,11 +542,14 @@ class ModifiedNFW(object): Parameters: log_Mhalo: float, optional - log10 of the Halo mass (solar masses) + log10 of the Halo mass [dark matter + baryons] (solar masses) + log_MCGM: float, optional + log10 of the CGM mass (solar masses) + If provided, this sets f_hot [this is its only use] c: float, optional concentration of the halo f_hot: float, optional - Fraction of the baryons in this hot phase + Fraction of the total "expected" baryons in this hot (aka CGM) phase Will likely use this for all diffuse gas alpha: float, optional Parameter to modify NFW profile power-law @@ -572,9 +573,21 @@ class ModifiedNFW(object): """ - def __init__(self, log_Mhalo=12.2, c=7.67, f_hot=0.75, alpha=0., - y0=1., z=0., cosmo=cosmo, **kwargs): + def __init__(self, + log_Mhalo=12.2, + c=7.67, + f_hot:float=None, #0.75, + log_MCGM:float=None, #1e12, + alpha=0., y0=1., z=0., cosmo=cosmo,fb = cosmo.Ob0/cosmo.Om0, + **kwargs): # Init + if f_hot is None: + if log_MCGM is None: + f_hot = 0.75 + log_MCGM = np.log10(f_hot * 10**log_Mhalo*fb) + else: + f_hot = 10**log_MCGM / (10**log_Mhalo*fb) + # Param self.log_Mhalo = log_Mhalo self.M_halo = 10.**self.log_Mhalo * constants.M_sun.cgs @@ -585,22 +598,22 @@ def __init__(self, log_Mhalo=12.2, c=7.67, f_hot=0.75, alpha=0., self.f_hot = f_hot self.zero_inner_ne = 0. # kpc self.cosmo = cosmo + self.fb = fb # Init more - self.setup_param(cosmo=self.cosmo) + self.setup_param(cosmo=self.cosmo,) - def setup_param(self,cosmo): + def setup_param(self,cosmo,): """ Setup key parameters of the model + + Args: + cosmo: astropy cosmology + Cosmology of the universe """ # Cosmology - if cosmo is None: - self.rhoc = 9.2e-30 * units.g / units.cm**3 - self.fb = 0.16 # Baryon fraction - self.H0 = 70. *units.km/units.s/ units.Mpc - else: - self.rhoc = self.cosmo.critical_density(self.z) - self.fb = cosmo.Ob0/cosmo.Om0 - self.H0 = cosmo.H0 + self.rhoc = self.cosmo.critical_density(self.z) + self.H0 = cosmo.H0 + # Dark Matter self.q = self.cosmo.Ode0/(self.cosmo.Ode0+self.cosmo.Om0*(1+self.z)**3) #r200 = (((3*Mlow*constants.M_sun.cgs) / (4*np.pi*200*rhoc))**(1/3)).to('kpc') @@ -1088,7 +1101,9 @@ class M31(ModifiedNFW): Taking mass from van der Marel 2012 """ - def __init__(self, log_Mhalo=12.18, c=7.67, f_hot=0.75, alpha=2, y0=2, **kwargs): + def __init__(self, + log_Mhalo=12.18, c=7.67, + alpha=2, y0=2,f_hot=0.75, **kwargs): # Init ModifiedNFW ModifiedNFW.__init__(self, log_Mhalo=log_Mhalo, c=c, f_hot=f_hot, diff --git a/frb/tests/test_frbhalos.py b/frb/tests/test_frbhalos.py index fbebe59c..62e9e7e6 100644 --- a/frb/tests/test_frbhalos.py +++ b/frb/tests/test_frbhalos.py @@ -87,6 +87,11 @@ def test_modified_NFW(): ne = mNFW.ne(xyz) assert np.all(ne > nH) +def test_log_MGCM(): + # Init + mNFW = halos.ModifiedNFW(log_MCGM=11.0) + assert np.isclose(mNFW.f_hot, 0.3989835640019133) + def test_milky_way(): Galaxy = halos.MilkyWay() assert np.isclose(Galaxy.M_halo.to('M_sun').value, 1.51356125e+12)