Skip to content

Commit

Permalink
take in units on vel
Browse files Browse the repository at this point in the history
  • Loading branch information
@xzackli
xzackli committed Feb 16, 2025
1 parent 09fe69f commit 1c2b951
Show file tree
Hide file tree
Showing 2 changed files with 17 additions and 141 deletions.
155 changes: 14 additions & 141 deletions src/profiles_ksz.jl
View file
Original file line number Diff line number Diff line change
@@ -1,27 +1,4 @@

# import PhysicalConstants.CODATA2018 as constants
# using Unitful
# const M_sun = 1.98847e30u"kg"
# const P_e_factor = constants.σ_e / (constants.m_e * constants.c_0^2)
# const T_cmb = 2.725 * u"K"
# using Cosmology
# using QuadGK
# using DelimitedFiles
# using Interpolations
# using NPZ


# table = npzread("/fs/lustre/cita/zack/projects/SZpack.v1.1.1/szpack_interp_4d.npz")
# nu_vector = LinRange(log(5.680062373019096*1e9),log(852.0093559528645*1e9),500)
# temp_vector = LinRange(1.0e-3,75.0,100)
# vel_vector = LinRange(1.0e-9, 1.0e-1,100)
# mu_vector = LinRange(-1,1,50)
# szpack_interp_ksz = scale(Interpolations.interpolate(table, BSpline(Cubic(Line(OnGrid())))), (temp_vector), (vel_vector), (mu_vector), (nu_vector))

# table_T0 = npzread("/fs/lustre/cita/zack/projects/SZpack.v1.1.1/szpack_interp_4d_0.npz")
# szpack_interp_T0 = scale(Interpolations.interpolate(table_T0[1,:,:,:], BSpline(Cubic(Line(OnGrid())))), (vel_vector), (mu_vector), (nu_vector))


struct Battaglia16KinematicSZProfile{T,C} <: AbstractGNFW{T}
f_b::T # Omega_b / Omega_c = 0.0486 / 0.2589
cosmo::C
Expand All @@ -37,7 +14,7 @@ function Battaglia16KinematicSZProfile(; Omega_c::T=0.2589, Omega_b::T=0.0486, h
end


struct RKSZpackProfile{T,C,ANG,P1,P2,ITP1,ITP2,ITP3} <: AbstractGNFW{T}
struct RKSZpackProfile{T,C,P1,P2,ITP1,ITP2,ITP3} <: AbstractGNFW{T}
f_b::T # Omega_b / Omega_c = 0.0486 / 0.2589
cosmo::C
X::T
Expand All @@ -48,13 +25,14 @@ struct RKSZpackProfile{T,C,ANG,P1,P2,ITP1,ITP2,ITP3} <: AbstractGNFW{T}
szpack_interp_T0::ITP3
end


function RKSZpackProfile(model_y::P1, model_tau::P2, tsz_interp::ITP1, szpack_interp_ksz::ITP2, szpack_interp_T0::ITP3;
Omega_c::T=0.2589, Omega_b::T=0.0486, h::T=0.6774, angle=true, x::T=2.6408) where {T <: Real, P1, P2, ITP1, ITP2, ITP3}
OmegaM=Omega_b+Omega_c
f_b = Omega_b / OmegaM
cosmo = get_cosmology(T, h=h, Neff=3.046, OmegaM=OmegaM)
return RKSZpackProfile{T, typeof(cosmo), angle, P1, P2, ITP1, ITP2, ITP3}(
f_b, cosmo, x, model_y, model_tau, tsz_interp, szpack_interp_ksz, szpack_interp_T0)
@assert isangletypeparameter(model_tau)
return RKSZpackProfile(f_b, cosmo, x, model_y, model_tau, tsz_interp, szpack_interp_ksz, szpack_interp_T0)
end


Expand All @@ -69,23 +47,24 @@ function SZpack_ksz(𝕡, M_200, z, r, vel; τ=0.01, mu = 1.0, showT=true)
θ_e = (constants.k_B*T_e)/(constants.m_e*constants.c_0^2)

# use velocity magnitude to determine direction along line-of-sight
if vel < 0
if sign(vel) < 0
mu *= -1
end

t = ustrip(uconvert(u"keV",T_e * constants.k_B))
nu = log(ustrip(X_to_nu(X)))
vel = abs(ustrip(vel/uconvert(u"km/s",constants.c_0))) # need to take absolute value of v to make sure v is within bounds of interpolator

vel = abs(uconvert(NoUnits, vel/constants.c_0))
# need to take absolute value of v to make sure v is within bounds of interpolator

# Term 1
dI_1 = uconvert(u"kg*s^-2",(𝕡.szpack_interp_ksz(t, vel, mu, nu) * u"MJy/sr" -
𝕡.szpack_interp_T0(vel, mu, nu) * u"MJy/sr")/τ)
y = XGPaint.compton_y(𝕡.p_y, M_200, z, r) #0
𝕡.szpack_interp_T0(vel, mu, nu) * u"MJy/sr") / τ)
y = XGPaint.compton_y(𝕡.p_y, M_200, z, r)
I_1 = uconvert(u"kg*s^-2",y * (dI_1/(θ_e)))

# Term 2
dI_2 = uconvert(u"kg*s^-2", (𝕡.szpack_interp_T0(vel, mu, nu) * u"MJy/sr")/τ)
tau = XGPaint.tau(𝕡.p_tau, r, M_200, z) #0.01 #XGPaint.tau_ksz(𝕡, M_200, z, r)
tau = XGPaint.tau(𝕡.p_tau, r, M_200, z)
I_2 = uconvert(u"kg*s^-2", dI_2 * tau)

I = I_1 + I_2
Expand Down Expand Up @@ -120,47 +99,9 @@ function non_rel_ksz(𝕡, M_200, z, r, vel; mu = 1.0)
end


function profile_grid_ksz(𝕡::AbstractGNFW{T}; N_z=256, N_logM=256, N_logθ=512, z_min=1e-3, z_max=5.0,
logM_min=11, logM_max=15.7, logθ_min=-16.5, logθ_max=2.5, N_v=256, v_min=-1353.6917, v_max=1276.7216) where T #logM_max=15.7

logθs = LinRange(logθ_min, logθ_max, N_logθ)
redshifts = LinRange(z_min, z_max, N_z)
logMs = LinRange(logM_min, logM_max, N_logM)
velocities = LinRange(v_min, v_max, N_v)

return profile_grid_ksz(𝕡, logθs, redshifts, logMs, velocities)
end


function profile_grid_ksz(𝕡::AbstractGNFW{T}, logθs, redshifts, logMs, velocities) where T

N_logθ, N_z, N_logM, N_vels = length(logθs), length(redshifts), length(logMs), length(velocities)
A = zeros(T, (N_logθ, N_z, N_logM, N_vels))

Threads.@threads for im in 1:N_logM
println("Completed Halo Mass $im")
logM = logMs[im]
M = 10^(logM) * M_sun
for (iz, z) in enumerate(redshifts)
for iθ in 1:N_logθ
θ = exp(logθs[iθ])
for (iv, vel) in enumerate(velocities)
szp = SZpack_ksz(𝕡, M, z, θ, vel) # FULL RELATIVISTIC
#szp = non_rel_ksz(𝕡, M, z, θ, vel) # NON-RELATIVISTIC
#A[iθ, iz, im, iv] = max(zero(T), szp)
A[iθ, iz, im, iv] = szp
end
end
end
end

return logθs, redshifts, logMs, velocities, A
end


function profile_paint_ksz!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}}, p,
function profile_paint!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}}, p::RKSZpackProfile,
α₀, δ₀, psa::CarClenshawCurtisProfileWorkspace,
sitp, z, Ms, vel, θmax) where T
z, Ms, vel, θmax) where T
# get indices of the region to work on
i1, j1 = sky2pix(m, α₀ - θmax, δ₀ - θmax)
i2, j2 = sky2pix(m, α₀ + θmax, δ₀ + θmax)
Expand Down Expand Up @@ -188,36 +129,7 @@ function profile_paint_ksz!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}}, p,
end
end


"""Helper function to build a tSZ interpolator"""
function build_interpolator_ksz(model::AbstractGNFW; cache_file::String="",
N_logθ=512, pad=256, overwrite=true, verbose=true)

if overwrite || (isfile(cache_file) == false)
verbose && print("Building new interpolator from model.\n")
rft = RadialFourierTransform(n=N_logθ, pad=pad)
logθ_min, logθ_max = log(minimum(rft.r)), log(maximum(rft.r))
prof_logθs, prof_redshift, prof_logMs, prof_velocity, prof_y = profile_grid_ksz(model;
N_logθ=N_logθ, logθ_min=logθ_min, logθ_max=logθ_max)
if length(cache_file) > 0
verbose && print("Saving new interpolator to $(cache_file).\n")
save(cache_file, Dict("prof_logθs"=>prof_logθs,
"prof_redshift"=>prof_redshift, "prof_logMs"=>prof_logMs, "prof_velocity"=>prof_velocity, "prof_y"=>prof_y))
end
else
print("Found cached Battaglia profile model. Loading from disk.\n")
model_grid = load(cache_file)
prof_logθs, prof_redshift, prof_logMs, prof_velocity, prof_y = model_grid["prof_logθs"],
model_grid["prof_redshift"], model_grid["prof_logMs"], model_grid["prof_velocity"], model_grid["prof_y"]
end

itp = Interpolations.interpolate(log.(prof_y), BSpline(Cubic(Line(OnGrid()))))
sitp = scale(itp, prof_logθs, prof_redshift, prof_logMs, prof_velocity)
return sitp
end


function profile_paint_ksz!(m::HealpixMap{T, RingOrder}, p,
function profile_paint!(m::HealpixMap{T, RingOrder}, p::RKSZpackProfile,
α₀, δ₀, w::HealpixProfileWorkspace, z, Mh, vel, θmax) where T
ϕ₀ = α₀
θ₀ = T(π)/2 - δ₀
Expand All @@ -239,42 +151,3 @@ function profile_paint_ksz!(m::HealpixMap{T, RingOrder}, p,
end
end


# for rectangular pixelizations

# multi-halo painting utilities
function paint_ksz!(m, p::XGPaint.AbstractProfile, psa, sitp,
masses::AV, redshifts::AV, αs::AV, δs::AV, velocities::AV, irange::AbstractUnitRange) where AV
for i in irange
α₀ = αs[i]
δ₀ = δs[i]
mh = masses[i]
z = redshifts[i]
v = velocities[i]
θmax_ = θmax(p, mh * XGPaint.M_sun, z)
profile_paint_ksz!(m, p, α₀, δ₀, psa, sitp, z, mh, v, θmax_)
end
end


function paint_ksz!(m, p::XGPaint.AbstractProfile, psa, sitp, masses::AV,
redshifts::AV, αs::AV, δs::AV, velocities::AV) where AV
fill!(m, 0)

N_sources = length(masses)
chunksize = ceil(Int, N_sources / (2Threads.nthreads()))
chunks = chunk(N_sources, chunksize);

Threads.@threads for i in 1:Threads.nthreads()
chunk_i = 2i
i1, i2 = chunks[chunk_i]
paint_ksz!(m, p, psa, sitp, masses, redshifts, αs, δs,velocities, i1:i2)
end

Threads.@threads for i in 1:Threads.nthreads()
chunk_i = 2i - 1
i1, i2 = chunks[chunk_i]
paint_ksz!(m, p, psa, sitp, masses, redshifts, αs, δs, velocities, i1:i2)
end
end

3 changes: 3 additions & 0 deletions src/tau.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,9 @@ struct BattagliaTauProfile{T,C,ANG} <: AbstractGNFW{T}
cosmo::C
end

# check if the tau profile is in terms of angle or distance; usually should be in angle
isangletypeparameter(::BattagliaTauProfile{T,C,true}) where {T,C} = true
isangletypeparameter(::BattagliaTauProfile{T,C,false}) where {T,C} = false

function BattagliaTauProfile(; Omega_c::T=0.2589, Omega_b::T=0.0486, h::T=0.6774, angle=true) where {T <: Real}
OmegaM=Omega_b+Omega_c
Expand Down

0 comments on commit 1c2b951

Please sign in to comment.