new value type for angle

src/tau.jl

@@ -1,14 +1,19 @@
 
-struct BattagliaTauProfile{T,C} <: AbstractGNFW{T}
+
+# ANG is a Value type, Val(true) if we want to use angular units, and 
+# Val(false) if we want to use physical units
+
+struct BattagliaTauProfile{T,C,ANG} <: AbstractGNFW{T}
     f_b::T  # Omega_b / Omega_c = 0.0486 / 0.2589
     cosmo::C
 end
 
-function BattagliaTauProfile(; Omega_c::T=0.2589, Omega_b::T=0.0486, h::T=0.6774) where {T <: Real}
+
+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
     f_b = Omega_b / OmegaM
     cosmo = get_cosmology(T, h=h, Neff=3.046, OmegaM=OmegaM)
-    return BattagliaTauProfile(f_b, cosmo)
+    return BattagliaTauProfile{T, typeof(cosmo), angle}(f_b, cosmo)
 end
 
 function get_params(::BattagliaTauProfile{T}, M_200, z) where T
@@ -33,19 +38,30 @@ function r200c_comoving(p, m200c, z)
 end
 
 
+function object_size(𝕡::BattagliaTauProfile{T,C,true}, physical_size, z) where {T,C}
+    d_A = angular_diameter_dist(𝕡.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
+
+function object_size(𝕡::BattagliaTauProfile{T,C,false}, physical_size, z) where {T,C}
+    return physical_size
+end
+
+
 # returns a density, which we can check against Msun/Mpc² 
-function rho_2d(p::BattagliaTauProfile, R_phys_projected, m200c, z)
+function rho_2d(p::BattagliaTauProfile, r, m200c, z)
     par = get_params(p, m200c, 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...
+    r200c = R_Δ(p, m200c, z, 200)
+    X = r / object_size(p, r200c, z)  # either ang/ang or phys/phys
+    rho_crit = ρ_crit_comoving_h⁻²(p, z)  # need to sort this out, it's all in comoving...
     result = par.P₀ * XGPaint._nfw_profile_los_quadrature(X, par.xc, par.α, par.β, par.γ)
 
     return result * rho_crit * (r200c * (1+z))
 end
 
-function ne2d(p::BattagliaTauProfile, R_phys_projected, m200c, z)
+function ne2d(p::BattagliaTauProfile, r, m200c, z)
     me = constants.ElectronMass
     mH = constants.ProtonMass
     mHe = 4mH
@@ -54,12 +70,14 @@ function ne2d(p::BattagliaTauProfile, R_phys_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_phys_projected, m200c, z)  # (Msun/h) / (Mpc/h)^2
+    result = rho_2d(p, r, m200c, z)  # (Msun/h) / (Mpc/h)^2
     return result / factor
 end
 
-function tau(p, R_phys_projected, m200c, z)
-    return constants.ThomsonCrossSection * ne2d(p, R_phys_projected, m200c, z) 
+# r is either a physical or angular radius. unitful does not do little h, so physical radius 
+# is always in Mpc, and angular radius is always in radians.
+function tau(p, r, m200c, z)
+    return constants.ThomsonCrossSection * ne2d(p, r, m200c, z) 
 end
 
 function profile_grid(𝕡::BattagliaTauProfile{T}, logθs, redshifts, logMs) where T

test/runtests.jl

@@ -140,7 +140,7 @@ end
 @testset "tau_profile" begin
     # test the 3D tau profile
 
-    p = BattagliaTauProfile(Omega_c=0.267, Omega_b=0.0493,  h=0.6712)
+    p = BattagliaTauProfile(Omega_c=0.267, Omega_b=0.0493,  h=0.6712, angle=false)
 
     # fits are in Msun/h, will change later
     par = get_params(p, (1e14M_sun / 0.6712), 0.5)