Cap HMCode2020 feedback stellar fraction at baryon fraction

Project.toml

@@ -28,3 +28,8 @@ LoopVectorization = "0.12.173"
 OrdinaryDiffEqCore = "3.25.0"
 OrdinaryDiffEqTsit5 = "1.10.0"
 Polyester = "0.7.19"
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]

src/feedback.jl

@@ -12,3 +12,20 @@ function get_feedback_parameters(T_AGN::Float64)
     )
     return params
 end
+
+"""
+    feedback_stellar_fraction(feedback_params, z, Omega_b, Omega_m)
+
+HMCode2020 feedback stellar fraction, capped at the cosmic baryon fraction.
+This follows the original Fortran HMcode implementation, where `HMcode_sbar`
+applies `fix_maximum(..., Omega_b/Omega_m)`.
+"""
+@inline function feedback_stellar_fraction(
+    feedback_params::Dict{Symbol, Float64},
+    z::Real,
+    Omega_b::Real,
+    Omega_m::Real,
+)
+    fstar = feedback_params[:f0] * 10.0^(float(z) * feedback_params[:fz])
+    return min(fstar, float(Omega_b) / float(Omega_m))
+end

src/power_spectrum.jl

@@ -225,7 +225,11 @@ function hmcode_power_single!(Pk_out, k, zs, cosmo, ws; T_AGN=nothing, tweaks=tr
         amp_fb = ws.M ./ rhom; compute_weights_inplace!(ws.w1h_mat_fb, ws.M, ws.nu_mat, amp_fb, ws.gcol_buf[1], ws.wtcol_buf[1]); w1h_ptr = ws.w1h_mat_fb
     end
     if (T_AGN !== nothing) && (!tweaks)
-        for iz in 1:nz; z = zs[iz]; ws.Mb_vec[iz] = feedback_params[:Mb0] * 10.0^(z * feedback_params[:Mbz]); ws.fstar_vec[iz] = feedback_params[:f0] * 10.0^(z * feedback_params[:fz]); end
+        for iz in 1:nz
+            z = zs[iz]
+            ws.Mb_vec[iz] = feedback_params[:Mb0] * 10.0^(z * feedback_params[:Mbz])
+            ws.fstar_vec[iz] = feedback_stellar_fraction(feedback_params, z, Om_b, Om_m)
+        end
     end
     for iz in 1:nz
         z = zs[iz]; dc, Dv, B = hmpars.delta_c[iz], hmpars.Delta_v[iz], hmpars.B[iz]

test/runtests.jl

@@ -5,6 +5,7 @@ using Test
     include("test_linear.jl")
     include("test_profiles.jl")
     include("test_halomodel.jl")
+    include("test_feedback.jl")
     include("test_power_spectrum.jl")
     include("test_10_cosmologies.jl")
 end

test/test_feedback.jl

@@ -0,0 +1,18 @@
+using Test
+include("../src/HMcode.jl")
+using .HMcode
+
+@testset "Feedback" begin
+    feedback_params = HMcode.get_feedback_parameters(10.0^7.7)
+    Omega_b = 0.051
+    Omega_m = 0.312
+    baryon_fraction = Omega_b / Omega_m
+
+    uncapped_low_z = feedback_params[:f0] * 10.0^(1.0 * feedback_params[:fz])
+    @test uncapped_low_z < baryon_fraction
+    @test HMcode.feedback_stellar_fraction(feedback_params, 1.0, Omega_b, Omega_m) ≈ uncapped_low_z
+
+    uncapped_high_z = feedback_params[:f0] * 10.0^(4.0 * feedback_params[:fz])
+    @test uncapped_high_z > baryon_fraction
+    @test HMcode.feedback_stellar_fraction(feedback_params, 4.0, Omega_b, Omega_m) == baryon_fraction
+end