Sigmoid

src/XGPaint.jl

@@ -38,6 +38,7 @@ include("./model.jl")
 include("./profiles.jl")
 include("./profiles_y.jl")
 include("./profiles_tau.jl")
+include("./profiles_tau_sigmoid.jl")
 include("./profiles_rsz.jl")
 include("./profiles_szp.jl")
 include("./profiles_rksz.jl")
@@ -51,6 +52,7 @@ export Radio_Sehgal2009, CIB_Planck2013, CIB_Scarfy, CO_CROWNED, LRG_Yuan23
 export paint!, generate_sources, process_sources, profile_grid, profile_paint!
 export profileworkspace, paint_szp!, profile_grid_szp, profile_paint_szp!, paint_rsz!, profile_grid_rsz, profile_paint_rsz!
 export build_interpolator, Battaglia16ThermalSZProfile, RSZPerturbativeProfile, build_interpolator_szp, build_interpolator_rsz
+export SigmoidBattagliaTauProfile, SigmoidBattagliaTauProfilePhysical, SigmoidBreakModel, sigmoid_suppression
 export SZPackRSZProfile, nu_to_X, X_to_nu, BattagliaTauProfile, HealpixRingProfileWorkspace
 
 end # module

src/profiles_tau_sigmoid.jl

@@ -0,0 +1,155 @@
+"""Sigmoid-suppressed Battaglia tau (optical depth) profiles.
+These profiles include a sigmoid suppression factor that smoothly transitions
+from 0.99 at 4 R_200 to 0.01 at 6 R_200.
+"""
+
+# can be used with either angular (default) or physical units
+#abstract type AbstractBattagliaTauProfile{T} <: AbstractGNFW{T} end
+
+# -------------------------------------------------------------------------
+# Sigmoid-suppressed tau models
+# -------------------------------------------------------------------------
+
+struct SigmoidBattagliaTauProfile{T,C} <: AbstractBattagliaTauProfile{T}
+    f_b::T
+    cosmo::C
+end
+
+struct SigmoidBattagliaTauProfilePhysical{T,C} <: AbstractBattagliaTauProfile{T}
+    f_b::T
+    cosmo::C
+end
+
+function SigmoidBattagliaTauProfile(; Omega_c::T=0.2589, Omega_b::T=0.0486, h::T=0.6774) where {T <: Real}
+    OmegaM = Omega_b + Omega_c
+    f_b    = Omega_b / OmegaM
+    cosmo  = get_cosmology(T, h=h, Neff=3.046, OmegaM=OmegaM)
+    return SigmoidBattagliaTauProfile{T, typeof(cosmo)}(f_b, cosmo)
+end
+
+function SigmoidBattagliaTauProfilePhysical(; Omega_c::T=0.2589, Omega_b::T=0.0486, h::T=0.6774) where {T <: Real}
+    OmegaM = Omega_b + Omega_c
+    f_b    = Omega_b / OmegaM
+    cosmo  = get_cosmology(T, h=h, Neff=3.046, OmegaM=OmegaM)
+    return SigmoidBattagliaTauProfilePhysical{T, typeof(cosmo)}(f_b, cosmo)
+end
+
+# if angular, return the R200 size in radians
+function object_size(model::SigmoidBattagliaTauProfile{T,C}, physical_size, z) where {T,C}
+    d_A = angular_diameter_dist(model.cosmo, z)
+    phys_siz_unitless = T(ustrip(uconvert(unit(d_A), physical_size)))
+    d_A_unitless = T(ustrip(d_A))
+    return atan(phys_siz_unitless, d_A_unitless)
+end
+
+# if physical, return the R200 size in Mpc
+function object_size(::SigmoidBattagliaTauProfilePhysical{T,C}, physical_size, z) where {T,C}
+    return physical_size
+end
+
+
+
+# -------------------------------------------------------------------------
+# Sigmoid definition
+# -------------------------------------------------------------------------
+"""
+    sigmoid_suppression(r)
+
+Sigmoid suppression factor that equals 0.99 at r=4 and 0.01 at r=6.
+Uses the form: f(r) = 1 / (1 + exp(k * (r - r0)))
+
+Where k and r0 are solved from the boundary conditions:
+- f(4) = 0.99 => 1/(1 + exp(k*(4-r0))) = 0.99
+- f(6) = 0.01 => 1/(1 + exp(k*(6-r0))) = 0.01
+
+This gives us:
+- exp(k*(4-r0)) = 1/99 ≈ 0.0101
+- exp(k*(6-r0)) = 99
+- Taking the ratio: exp(2k) = 99^2 = 9801
+- So k = ln(9801)/2 ≈ 4.599
+- And r0 = 4 + ln(99)/k = 4 + ln(99)/4.599 ≈ 5
+"""
+function sigmoid_suppression(r::T) where {T <: Real}
+    k  = T(2.3)
+    r0 = T(5.0)
+    return 1 / (1 + exp(k * (r - r0)))
+end
+
+# -------------------------------------------------------------------------
+# Suppressed line-of-sight integration
+# -------------------------------------------------------------------------
+function _sigmoid_nfw_profile_los_quadrature(x, xc, α, β, γ; zmax=1e5, rtol=eps(), order=9)
+    x² = x^2
+    scale = 1e9
+    integral, err = quadgk(
+        y -> scale * sigmoid_suppression(√(y^2 + x²)) *
+             XGPaint.generalized_nfw(√(y^2 + x²), xc, α, β, γ),
+        0.0, zmax, rtol=rtol, order=order)
+    return 2integral / scale
+end
+
+# -------------------------------------------------------------------------
+# Density and number profiles with sigmoid suppression
+# -------------------------------------------------------------------------
+function rho_2d(model::SigmoidBattagliaTauProfile, r, m200c, z)
+    par     = get_params(model, m200c, z)
+    r200c   = R_Δ(model, m200c, z, 200)
+    X       = r / object_size(model, r200c, z)
+    rho_crit = ρ_crit_comoving_h⁻²(model, z)
+    result  = par.P₀ * _sigmoid_nfw_profile_los_quadrature(X, par.xc, par.α, par.β, par.γ)
+    return result * rho_crit * (r200c * (1 + z))
+end
+
+function rho_2d(model::SigmoidBattagliaTauProfilePhysical, r, m200c, z)
+    par     = get_params(model, m200c, z)
+    r200c   = R_Δ(model, m200c, z, 200)
+    X       = r / object_size(model, r200c, z)
+    rho_crit = ρ_crit_comoving_h⁻²(model, z)
+    result  = par.P₀ * _sigmoid_nfw_profile_los_quadrature(X, par.xc, par.α, par.β, par.γ)
+    return result * rho_crit * (r200c * (1 + z))
+end
+"""
+function ne2d(model::SigmoidBattagliaTauProfile, r, m200c, z)
+    me  = constants.ElectronMass
+    mH  = constants.ProtonMass
+    mHe = 4mH
+    xH  = 0.76
+    nH_ne  = 2xH / (xH + 1)
+    nHe_ne = (1 - xH)/(2 * (1 + xH))
+    factor = (me + nH_ne*mH + nHe_ne*mHe) / model.cosmo.h^2
+    result = rho_2d(model, r, m200c, z)
+    return result / factor
+end
+
+function ne2d(model::SigmoidBattagliaTauProfilePhysical, r, m200c, z)
+    me  = constants.ElectronMass
+    mH  = constants.ProtonMass
+    mHe = 4mH
+    xH  = 0.76
+    nH_ne  = 2xH / (xH + 1)
+    nHe_ne = (1 - xH)/(2 * (1 + xH))
+    factor = (me + nH_ne*mH + nHe_ne*mHe) / model.cosmo.h^2
+    result = rho_2d(model, r, m200c, z)
+    return result / factor
+end
+"""
+# -------------------------------------------------------------------------
+# Compute τ (optical depth)
+# -------------------------------------------------------------------------
+"""
+function compute_tau(model::SigmoidBattagliaTauProfile, r, m200c, z)
+    return constants.ThomsonCrossSection * ne2d(model, r, m200c, z) + 0
+end
+
+function compute_tau(model::SigmoidBattagliaTauProfilePhysical, r, m200c, z)
+    return constants.ThomsonCrossSection * ne2d(model, r, m200c, z) + 0
+end
+"""
+# -------------------------------------------------------------------------
+# Direct call overloads
+# -------------------------------------------------------------------------
+#(model::SigmoidBattagliaTauProfile)(r, m200c, z) =
+#    compute_tau(model, r, m200c * M_sun, z)
+
+#(model::SigmoidBattagliaTauProfilePhysical)(r, m200c, z) =
+#    compute_tau(model, r, m200c * M_sun, z)

test/test_sigmoid.jl

@@ -0,0 +1,94 @@
+#!/usr/bin/env julia
+
+# Test script to evaluate and plot signmoid τ (optical depth) profiles
+
+using XGPaint
+using Plots
+using Printf
+
+# ------------------------------------------------------------
+# Test parameters - typical cluster
+# ------------------------------------------------------------
+M_test = 1e14   # solar masses (unitless; multiplied by M_sun internally)
+z_test = 0.3
+
+println("Creating τ-profiles for test cluster:")
+println("  Mass: $(M_test) M_sun")
+println("  Redshift: $(z_test)")
+
+# ------------------------------------------------------------
+# Create both standard and sigmoid τ models
+# ------------------------------------------------------------
+tau_model = BattagliaTauProfile()
+sigmoid_tau_model = SigmoidBattagliaTauProfile()
+
+println("Models created successfully!")
+
+# ------------------------------------------------------------
+# Compute R_200 and angular size
+# ------------------------------------------------------------
+R_200_test = XGPaint.R_Δ(tau_model, M_test * XGPaint.M_sun, z_test, 200)
+angular_size_test = XGPaint.object_size(tau_model, R_200_test, z_test)
+
+println("Cluster properties:")
+println("  R_200 = $(R_200_test / 1000) kpc")  # approximate conversion
+println("  Angular size of R_200 = $(angular_size_test * 180/π * 3600) arcsec")
+
+# ------------------------------------------------------------
+# Define radial range: 0.01–10 R200 (in angular units)
+# ------------------------------------------------------------
+r_over_R200 = 10 .^ range(log10(0.01), log10(10), length=100)
+radii = r_over_R200 .* angular_size_test   # angular radii in radians
+
+# ------------------------------------------------------------
+# Evaluate τ(r)
+# ------------------------------------------------------------
+println("Evaluating standard Battaglia τ profile...")
+tau_standard = [tau_model(r, M_test, z_test) for r in radii]
+
+println("Evaluating sigmoid-suppressed τ profile...")
+tau_sigmoid = [sigmoid_tau_model(r, M_test, z_test) for r in radii]
+
+println("Profile evaluation completed!")
+
+# ------------------------------------------------------------
+# Plot with Plots.jl (GR backend)
+# ------------------------------------------------------------
+gr()
+plot(r_over_R200, tau_standard;
+     lw=2, xscale=:log10, yscale=:log10,
+     label="Standard Battaglia τ", color=:blue)
+plot!(r_over_R200, tau_sigmoid;
+      lw=2, label="Sigmoid-Suppressed τ", color=:red,
+      xlabel="r / R₂₀₀", ylabel="Optical Depth τ(r)",
+      title="Optical Depth Profiles: Standard vs Sigmoid-Suppressed\nM = 10¹⁴ M☉, z = 0.3")
+
+# Add dashed lines for the sigmoid region
+vline!([4.0, 6.0], line=(:dash, :gray), label=["r=4R₂₀₀" "r=6R₂₀₀"])
+
+# Save the plot
+savefig("tau_profiles_comparison.png")
+println("Plot saved as 'tau_profiles_comparison.png'")
+
+# ------------------------------------------------------------
+# Print summary table
+# ------------------------------------------------------------
+println("\nProfile comparison at key radii:")
+println("r/R₂₀₀    Standard τ       Sigmoid τ        Suppression")
+println("--------------------------------------------------------")
+for r_ratio in [0.1, 1.0, 4.0, 5.0, 6.0, 8.0]
+    idx = argmin(abs.(r_over_R200 .- r_ratio))
+    t_std = tau_standard[idx]
+    t_sig = tau_sigmoid[idx]
+    suppression = t_sig / t_std
+    @printf("%.1f       %.2e      %.2e      %.3f\n", r_ratio, t_std, t_sig, suppression)
+end
+
+# ------------------------------------------------------------
+# Check sigmoid suppression function directly
+# ------------------------------------------------------------
+println("\nSigmoid suppression function values:")
+for r_ratio in [3.0, 4.0, 5.0, 6.0, 7.0]
+    supp = sigmoid_suppression(r_ratio)
+    @printf("sigmoid_suppression(%.1f) = %.4f\n", r_ratio, supp)
+end