working now

Author: alexander-mead

README.md

@@ -33,6 +33,7 @@ import hmcode
 # Ranges
 k = np.logspace(-3, 1, 100) # Wavenumbers [h/Mpc]
 zs = [3., 2., 1., 0.5, 0.]  # Redshifts
+T_AGN = 10**7.8             # Feedback temperature [K]
 
 # Run CAMB
 parameters = camb.CAMBparams(WantCls=False)
@@ -41,7 +42,7 @@ parameters.set_matter_power(redshifts=zs, kmax=100.) # kmax should be much large
 results = camb.get_results(parameters)
 
 # HMcode
-Pk = hmcode.power(k, zs, results)
+Pk = hmcode.power(k, zs, results, T_AGN)
 ```
 
 ## Note
@@ -80,7 +81,7 @@ I think any residual differences between codes must therefore stem from:
 But I didn't have time to investigate these differences more thoroughly. Note that there are accuracy parameters in `CAMB-HMcode` fixed at the $10^{-4}$ level, so you would never expect better than 0.01% agreement. Given that `HMcode` is only accurate at the ~2.5% level compared to simulated power spectra, the level of agreement between the codes seems okay to me, with the caveats above regarding very massive neutrinos.
 
 While writing this code I had a few ideas for future improvements:
-- Add the `HMcode-2020` baryon-feedback model; this would not be too hard for the enthusiastic student/postdoc.
+- ~~Add the `HMcode-2020` baryon-feedback model; this would not be too hard for the enthusiastic student/postdoc.~~ (Thanks Laila Linke!)
 - The predictions are a bit sensitive to the smoothing $\sigma$ used for the dewiggling. This should probably be a fitted parameter.
 - It's annoying having to calculate linear growth functions (all, LCDM), especially since the linear growth doesn't really exist. One should probably use the $P(k)$ amplitude evolution over time at some cleverly chosen scale instead, or instead the evolution of $\sigma(R)$ over time at some pertinent $R$ value. Note that the growth factors are *only* used to calculate the [Dolag et al. (2004)](https://arxiv.org/abs/astro-ph/0309771) correction and [Mead (2017)](https://arxiv.org/abs/1606.05345) $\delta_\mathrm{c}$, $\Delta_\mathrm{v}$.
 - I never liked the halo bloating parameter, it's hard to understand the effect of modifying halo profiles in Fourier space. Someone should get rid of this (maybe modify the halo mass function instead?).

hmcode/__init__.py

@@ -1,2 +1 @@
-from .hmcode import power
-from .hmcode import _get_feedback_suppression as get_feedback_suppression
+from .hmcode import power

hmcode/hmcode.py

@@ -157,7 +157,7 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, T_AGN=None,
         else: # Use vanilla-ish halo model if no HMcode tweaks used
             B, eta = 4., 0.
 
-        if T_AGN:
+        if T_AGN and not tweaks:
             B = feedback_params['B0']*np.power(10, z*feedback_params['Bz'])
             Mb = feedback_params['Mb0']*np.power(10, z*feedback_params['Mbz'])
             fstar = feedback_params['f0']*np.power(10, z*feedback_params['fz'])
@@ -171,7 +171,7 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, T_AGN=None,
                 print('Halo bloating; eta: {:.4}'.format(eta))
                 print('Minimum halo concentration; B: {:.4}'.format(B))#, B/4.)
                 print('Transition smoothing; alpha: {:.4}'.format(alpha))
-            if T_AGN:
+            if T_AGN and not tweaks:
                 print('Gas-loss halo-mass parameter Mb: {:1e}'.format(Mb))
                 print('Effective halo stellar-mass fraction f*: {:.4f}'.format(f))
             print()
@@ -186,11 +186,11 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, T_AGN=None,
         rv = hmod.virial_radius(M) 
         Uk = np.ones((len(k), len(M)))
         for iM, (_rv, _c, _nu) in enumerate(zip(rv, c, nu)): # TODO: Remove loop for speed?
-            if T_AGN:
+            if T_AGN and not tweaks:
                 Uk[:, iM] = _win_NFW_baryons(k*(_nu**eta), _rv, _c, M[iM], Mb, fstar, Om_m, Om_c, Om_b)[:, 0]
             else:
                 Uk[:, iM] = _win_NFW(k*(_nu**eta), _rv, _c)[:, 0]
-        if T_AGN: # NOTE: No Factor of 1-f_nu for baryonic, because this effect is already included!
+        if T_AGN and not tweaks: # NOTE: No Factor of 1-f_nu for baryonic, because this effect is already included!
             profile = halo.profile.Fourier(k, M, Uk, amplitude=M/hmod.rhom, mass_tracer=True) 
         else: # NOTE: Factor of 1-f_nu in profile amplitude
             profile = halo.profile.Fourier(k, M, Uk, amplitude=M*(1.-f_nu)/hmod.rhom, mass_tracer=True) 
@@ -209,7 +209,7 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, T_AGN=None,
             Pk_hm = Pk_lin+_Pk_1h['m-m']
         Pk_HMcode[iz, :] = Pk_hm
 
-    if tweaks and T_AGN:
+    if T_AGN and tweaks:
         suppression = _get_feedback_suppression(k, zs, CAMB_results, T_AGN, Mmin=Mmin, Mmax=Mmax, nM=nM, verbose=False)
         Pk_HMcode *= suppression