switch to taking in physical radius

WebSky-CITA · Nov 8, 2024 · 42db07c · 42db07c
1 parent 5029497
commit 42db07c
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 33 deletions.
diff --git a/docs/src/ksz.md b/docs/src/ksz.md
@@ -48,8 +48,9 @@ p = BattagliaTauProfile(Omega_c=0.267, Omega_b=0.0493,  h=0.6712)
 We can now compute the integrated electron density, 
 
 ```@example ksz
+using Unitful, UnitfulAstro
 zz = 2.5
-Mnew = 93218298413772.23 * (M_sun / p.cosmo.h)  # Mcrit
+Mnew = 93218298413772.23 * (XGPaint.M_sun / p.cosmo.h)  # Mcrit
 tc = ne2d(p, 3u"Mpc" / p.cosmo.h, Mnew, zz) / (p.cosmo.h / 1u"Mpc")^2 + 0
 @test abs(1 - tc / 1.4217533501414173e+69) < 1e-3
 ```

diff --git a/src/tau.jl b/src/tau.jl
@@ -22,25 +22,6 @@ function get_params(::BattagliaTauProfile{T}, M_200, z) where T
     return (xc=T(xc), α=T(α), β=T(β), γ=T(γ), P₀=T(P₀))
 end
 
-function profile_grid(𝕡::BattagliaTauProfile{T}, logθs, redshifts, logMs) where T
-
-    N_logθ, N_z, N_logM = length(logθs), length(redshifts), length(logMs)
-    A = zeros(T, (N_logθ, N_z, N_logM))
-
-    Threads.@threads :static for im in 1:N_logM
-        logM = logMs[im]
-        M = 10^(logM) * M_sun
-        for (iz, z) in enumerate(redshifts)
-            for iθ in 1:N_logθ
-                θ = exp(logθs[iθ])
-                y = compton_y(𝕡, M, z, θ)
-                A[iθ, iz, im] = max(zero(T), y)
-            end
-        end
-    end
-
-    return logθs, redshifts, logMs, A
-end
 
 function ρ_crit_comoving_h⁻²(p, z)
     return  (ρ_crit(p, z) ) / (1+z)^3 / p.cosmo.h^2
@@ -53,17 +34,18 @@ end
 
 
 # returns a density, which we can check against Msun/Mpc² 
-function rho_2d(p::BattagliaTauProfile, R_comoving_projected, m200c, z)
+function rho_2d(p::BattagliaTauProfile, R_phys_projected, m200c, z)
     par = get_params(p, m200c, z)
-    r200c = r200c_comoving(p, m200c, z)
-    X = R_comoving_projected / r200c
-    rho_crit = ρ_crit_comoving_h⁻²(p, z)
+    # r200c = r200c_comoving(p, m200c, z)
+    r200c = R_Δ(p, m200c, z, 200) 
+    X = R_phys_projected/ (r200c) # X is calculated in phys / phys
+    rho_crit = ρ_crit_comoving_h⁻²(p, z)  # need to sort this out; calc is in comoving...
     result = par.P₀ * XGPaint._nfw_profile_los_quadrature(X, par.xc, par.α, par.β, par.γ)
 
-    return result * rho_crit * r200c
+    return result * rho_crit * (r200c * (1+z))
 end
 
-function ne2d(p::BattagliaTauProfile, R_comoving_projected, m200c, z)
+function ne2d(p::BattagliaTauProfile, R_phys_projected, m200c, z)
     me = constants.ElectronMass
     mH = constants.ProtonMass
     mHe = 4mH
@@ -72,12 +54,12 @@ function ne2d(p::BattagliaTauProfile, R_comoving_projected, m200c, z)
     nHe_ne = (1 - xH)/(2 * (1 + xH))
     factor = (me + nH_ne*mH + nHe_ne*mHe) / p.cosmo.h^2
 
-    result = rho_2d(p, R_comoving_projected, m200c, z)  # (Msun/h) / (Mpc/h)^2
+    result = rho_2d(p, R_phys_projected, m200c, z)  # (Msun/h) / (Mpc/h)^2
     return result / factor
 end
 
-function tau(p, R_comoving_projected, m200c, z)
-    return constants.ThomsonCrossSection * ne2d(p, R_comoving_projected, m200c, z) 
+function tau(p, R_phys_projected, m200c, z)
+    return constants.ThomsonCrossSection * ne2d(p, R_phys_projected, m200c, z) 
 end
 
 function profile_grid(𝕡::BattagliaTauProfile{T}, logθs, redshifts, logMs) where T

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -159,21 +159,25 @@ end
     zz = 0.5
 
     Mnew = 86349927525539.45 * (M_sun / p.cosmo.h)
-    ta = rho_2d(p, 1u"Mpc" / p.cosmo.h, Mnew, zz) / (M_sun / 1u"Mpc^2") * p.cosmo.h + 0
+    ta = rho_2d(p, 1u"Mpc" / (1+zz) / p.cosmo.h, Mnew, zz) / (M_sun / 1u"Mpc^2") * p.cosmo.h + 0
     @test abs(1 - ta / 10.149657232478662e12) < 1e-4
 
-    tb = rho_2d(p, 2u"Mpc" / p.cosmo.h, Mnew, zz) / (M_sun / 1u"Mpc^2") * p.cosmo.h + 0
+    tb = rho_2d(p, 2u"Mpc" / (1+zz) / p.cosmo.h, Mnew, zz) / (M_sun / 1u"Mpc^2") * p.cosmo.h + 0
     @test abs(1 - tb / 2.446742995577804e12) < 1e-4
 
     zz = 2.5
     Mnew = 93218298413772.23 * (M_sun / p.cosmo.h)  # Mcrit
-    tc = rho_2d(p, 3u"Mpc" / p.cosmo.h, Mnew, zz) / (M_sun / 1u"Mpc^2") * p.cosmo.h + 0
+    tc = rho_2d(p, 3u"Mpc" / (1+zz) / p.cosmo.h, Mnew, zz) / (M_sun / 1u"Mpc^2") * p.cosmo.h + 0
     @test abs(1 - tc / 0.9123565449059163e12) < 1e-4
 
     zz = 2.5
     Mnew = 93218298413772.23 * (M_sun / p.cosmo.h)  # Mcrit
-    tc = ne2d(p, 3u"Mpc" / p.cosmo.h, Mnew, zz) / (p.cosmo.h / 1u"Mpc")^2 + 0
+    tc = ne2d(p, 3u"Mpc" / (1+zz) / p.cosmo.h, Mnew, zz) / (p.cosmo.h / 1u"Mpc")^2 + 0
     @test abs(1 - tc / 1.4217533501414173e+69) < 1e-3
 
+    zz = 2.5
+    Mnew = 93218298413772.23 * (XGPaint.M_sun / p.cosmo.h)  # Mcrit
+    tc = XGPaint.tau(p, 3u"Mpc" / (1+zz) / p.cosmo.h, Mnew, zz) + 0
+    @test abs(1 - tc / 4.475127577749756e-05) < 1e-3
 
 end