diff --git a/autogalaxy/config/priors/mass/dark/yang24.yaml b/autogalaxy/config/priors/mass/dark/yang24.yaml new file mode 100644 index 00000000..f2ef8589 --- /dev/null +++ b/autogalaxy/config/priors/mass/dark/yang24.yaml @@ -0,0 +1,51 @@ +YangSIDMSph: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + kappa_s: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.2 + limits: + lower: 0.0 + upper: inf + scale_radius: + type: Uniform + lower_limit: 0.0 + upper_limit: 30.0 + width_modifier: + type: Relative + value: 0.2 + limits: + lower: 0.0 + upper: inf + tau: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.2 + limits: + lower: 0.0 + upper: 1.0 diff --git a/autogalaxy/config/priors/mass/dark/yang24_mcr.yaml b/autogalaxy/config/priors/mass/dark/yang24_mcr.yaml new file mode 100644 index 00000000..df10d937 --- /dev/null +++ b/autogalaxy/config/priors/mass/dark/yang24_mcr.yaml @@ -0,0 +1,91 @@ +YangSIDMMCRLudlowSph: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + mass_at_200: + type: LogUniform + lower_limit: 100000000.0 + upper_limit: 1000000000000000.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf + sigma_over_m: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf + velocity_exponent: + type: Uniform + lower_limit: 0.0 + upper_limit: 4.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: 8.0 + velocity_ref: + type: Uniform + lower_limit: 1.0 + upper_limit: 100.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf + t_age: + type: Uniform + lower_limit: 0.0 + upper_limit: 13.8 + width_modifier: + type: Relative + value: 0.2 + limits: + lower: 0.0 + upper: inf + redshift_object: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf + redshift_source: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py index 19b4cd8f..e724bef0 100644 --- a/autogalaxy/profiles/mass/__init__.py +++ b/autogalaxy/profiles/mass/__init__.py @@ -48,6 +48,8 @@ cNFWMCRScatterLudlowSph, KaplinghatCoredNFWSph, KaplinghatCoredNFWMCRLudlowSph, + YangSIDMSph, + YangSIDMMCRLudlowSph, ) from .stellar import ( Gaussian, diff --git a/autogalaxy/profiles/mass/dark/__init__.py b/autogalaxy/profiles/mass/dark/__init__.py index dac908e7..63afe797 100644 --- a/autogalaxy/profiles/mass/dark/__init__.py +++ b/autogalaxy/profiles/mass/dark/__init__.py @@ -14,3 +14,5 @@ from .cnfw_mcr_scatter import cNFWMCRScatterLudlow, cNFWMCRScatterLudlowSph from .kaplinghat import KaplinghatCoredNFWSph from .kaplinghat_mcr import KaplinghatCoredNFWMCRLudlowSph +from .yang24 import YangSIDMSph +from .yang24_mcr import YangSIDMMCRLudlowSph diff --git a/autogalaxy/profiles/mass/dark/yang24.py b/autogalaxy/profiles/mass/dark/yang24.py new file mode 100644 index 00000000..4e06d385 --- /dev/null +++ b/autogalaxy/profiles/mass/dark/yang24.py @@ -0,0 +1,296 @@ +from typing import Tuple + +import numpy as np +from scipy.integrate import quad + +import autoarray as aa + +from autogalaxy.profiles.mass.dark.abstract import DarkProfile +from autogalaxy.profiles.mass.dark.nfw import NFWSph +from autogalaxy.profiles.mass.dark.kaplinghat import ( + _nfw_radial_deflection_from, + _trapezoid_from, +) +from autogalaxy.profiles.mass.abstract.abstract import MassProfile + + +def _yang24_parameter_ratios_from(tau): + """ + The Yang et al. (2024) fitting functions for how the cored-NFW parameters + evolve with the dimensionless gravothermal time ``tau = t / t_c`` + (arXiv:2305.16176, Section 2). + + Returns the ratios (rho_s / rho_s0, r_s / r_s0, r_c / r_s0) relative to the + halo's CDM baseline NFW parameters (rho_s0, r_s0). All three fits are exact + at ``tau = 0`` (ratios 1, 1, 0), recovering the NFW profile. + """ + tau = float(tau) + + log_term = np.log(tau + 0.001) / np.log(0.001) + + rho_s_ratio = ( + 2.033 + + 0.7381 * tau + + 7.264 * tau**5 + - 12.73 * tau**7 + + 9.915 * tau**9 + + (1.0 - 2.033) * log_term + ) + r_s_ratio = ( + 0.7178 + - 0.1026 * tau + + 0.2474 * tau**2 + - 0.4079 * tau**3 + + (1.0 - 0.7178) * log_term + ) + r_c_ratio = ( + 2.555 * np.sqrt(tau) + - 3.632 * tau + + 2.131 * tau**2 + - 1.415 * tau**3 + + 0.4683 * tau**4 + ) + + return rho_s_ratio, r_s_ratio, max(r_c_ratio, 0.0) + + +def _yang24_density_3d_from(radii, rho_s, r_s, r_c, xp): + """ + The beta = 4 cored-NFW density of Yang et al. (2024): + + rho(r) = rho_s / [ ((r^4 + r_c^4)^(1/4) / r_s) (1 + r / r_s)^2 ] + + Closed form in any array backend; reduces exactly to NFW as r_c -> 0. + """ + r = xp.maximum(radii, 1.0e-12) + r_tilde = (r**4 + r_c**4) ** 0.25 + return rho_s / ((r_tilde / r_s) * (1.0 + r / r_s) ** 2) + + +class YangSIDMSph(MassProfile, DarkProfile): + r""" + Spherical SIDM halo following the parametric gravothermal-evolution model of + Yang, Nadler, Yu & Zhong (2024), arXiv:2305.16176. + + The density is the beta = 4 cored-NFW form + + .. math:: + \rho(r) = \frac{\rho_s}{\left[(r^4 + r_c^4)^{1/4} / r_s\right] + \left(1 + r / r_s\right)^2} + + whose parameters :math:`(\rho_s, r_s, r_c)` evolve along the universal + gravothermal sequence as a function of the dimensionless time + :math:`\tau = t / t_c`, via the Yang et al. (2024) fitting functions. At + ``tau = 0`` the profile is exactly the NFW profile of the halo's CDM + baseline ``(kappa_s, scale_radius)``; as ``tau`` grows the halo first forms + a core (core expansion) and then contracts towards core collapse. + + The fits are calibrated for ``tau`` in [0, 1]; input values are truncated + into that range, following the truncation applied in Yang et al. (2024) + for deep collapse. + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + kappa_s + The overall normalization of the CDM baseline NFW halo + (rho_s0 * scale_radius / critical surface density). + scale_radius + The arc-second scale radius r_s0 of the CDM baseline NFW halo. + tau + The dimensionless gravothermal evolution time t / t_c, where t_c is the + halo's core-collapse timescale. Use `YangSIDMMCRLudlowSph` to compute + ``tau`` from physical inputs (cross section, halo age and mass). + """ + + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + kappa_s: float = 0.05, + scale_radius: float = 1.0, + tau: float = 0.0, + ): + super().__init__(centre=centre, ell_comps=(0.0, 0.0)) + + self.kappa_s = kappa_s + self.scale_radius = scale_radius + self.tau = min(max(float(tau), 0.0), 1.0) + + rho_s_ratio, r_s_ratio, r_c_ratio = _yang24_parameter_ratios_from( + tau=self.tau + ) + + self.rho_s_evolved = rho_s_ratio * kappa_s / scale_radius + self.scale_radius_evolved = r_s_ratio * scale_radius + self.core_radius_evolved = r_c_ratio * scale_radius + + self._nfw = NFWSph( + centre=(0.0, 0.0), + kappa_s=kappa_s, + scale_radius=scale_radius, + ) + + def _density_3d_from_radius(self, radii): + radii = np.asarray(radii, dtype=float) + return _yang24_density_3d_from( + radii=radii, + rho_s=self.rho_s_evolved, + r_s=self.scale_radius_evolved, + r_c=self.core_radius_evolved, + xp=np, + ) + + def density_3d_func(self, r, xp=np): + radii = r.array if hasattr(r, "array") else r + if xp is not np: + radii = np.asarray(radii) + return self._density_3d_from_radius(radii) + + def convergence_func(self, grid_radius, xp=np): + radii = ( + grid_radius.array + if hasattr(grid_radius, "array") + else np.asarray(grid_radius) + ) + scalar_input = np.ndim(radii) == 0 + radii = np.atleast_1d(np.asarray(radii, dtype=float)) + + if self.tau <= 0.0: + values = self._nfw.convergence_func(aa.ArrayIrregular(radii), xp=np) + return values[0] if scalar_input else values + + z_max = max( + 500.0 * self.scale_radius_evolved, 50.0 * self.core_radius_evolved + ) + + def convergence_at_radius(radius): + radius = float(max(radius, 1.0e-8)) + integral = quad( + lambda z: self._density_3d_from_radius(np.sqrt(radius**2 + z**2)), + 0.0, + z_max, + epsrel=1.0e-5, + limit=100, + )[0] + return 2.0 * integral + + convergence = np.array([convergence_at_radius(radius) for radius in radii]) + return convergence[0] if scalar_input else convergence + + @aa.over_sample + @aa.decorators.to_array + @aa.decorators.transform + def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): + radii = self.radial_grid_from(grid=grid, xp=np, **kwargs) + return self.convergence_func(grid_radius=radii, xp=np) + + def radial_deflection_from_radius(self, radius): + radius = float(radius) + + if radius <= 1.0e-8: + return 0.0 + + if self.tau <= 0.0: + return float( + np.sqrt( + np.sum( + self._nfw.deflections_yx_2d_from( + grid=aa.Grid2DIrregular([[radius, 0.0]]) + ).array[0] + ** 2.0 + ) + ) + ) + + mass_2d = quad( + lambda r: self.convergence_func(aa.ArrayIrregular([r]))[0] * r, + 0.0, + radius, + epsrel=1.0e-4, + limit=100, + )[0] + return 2.0 * mass_2d / radius + + @staticmethod + def radial_deflection_from(r, params, xp): + rho_s = params[0] + r_s = params[1] + r_c = params[2] + + r = xp.asarray(r) + r_safe = xp.maximum(r, 1.0e-8) + + z_max = xp.maximum(500.0 * r_s, 50.0 * r_c) + z_unit = xp.linspace(1.0e-5, 1.0, 160) + z = z_max * z_unit**3 + u = xp.linspace(0.0, 1.0, 64) + + projected_radii = xp.maximum(r_safe[:, None] * u[None, :], 1.0e-6) + three_d_radii = xp.sqrt( + projected_radii[:, :, None] ** 2 + z[None, None, :] ** 2 + ) + + density = _yang24_density_3d_from( + radii=three_d_radii, + rho_s=rho_s, + r_s=r_s, + r_c=r_c, + xp=xp, + ) + convergence = 2.0 * _trapezoid_from(density, x=z, axis=-1, xp=xp) + + mass_integral = r_safe**2 * _trapezoid_from( + convergence * u[None, :], x=u, axis=-1, xp=xp + ) + numerical = 2.0 * mass_integral / r_safe + analytic_nfw = _nfw_radial_deflection_from( + r=r_safe, + kappa_s=rho_s * r_s, + scale_radius=r_s, + xp=xp, + ) + + radial_deflection = xp.where( + r_c > 0.0, + numerical, + analytic_nfw, + ) + + return xp.where(r > 1.0e-8, radial_deflection, 0.0) + + @aa.decorators.to_vector_yx + @aa.decorators.transform + def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): + theta = self.radial_grid_from(grid=grid, xp=np, **kwargs).array + deflection_r = np.array( + [self.radial_deflection_from_radius(radius) for radius in theta] + ) + + return self._cartesian_grid_via_radial_from( + grid=grid, + radius=deflection_r, + xp=np, + **kwargs, + ) + + @aa.over_sample + @aa.decorators.to_array + @aa.decorators.transform + def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): + theta = self.radial_grid_from(grid=grid, xp=np, **kwargs).array + + potential = np.array( + [ + quad( + self.radial_deflection_from_radius, + 0.0, + max(float(radius), 1.0e-8), + epsrel=1.0e-4, + limit=100, + )[0] + for radius in theta + ] + ) + + return potential diff --git a/autogalaxy/profiles/mass/dark/yang24_mcr.py b/autogalaxy/profiles/mass/dark/yang24_mcr.py new file mode 100644 index 00000000..59bbed7f --- /dev/null +++ b/autogalaxy/profiles/mass/dark/yang24_mcr.py @@ -0,0 +1,228 @@ +import math +from typing import Tuple + +import numpy as np + +from autogalaxy.profiles.mass.dark import mcr_util +from autogalaxy.profiles.mass.dark.yang24 import YangSIDMSph + + +MSUN_G = 1.98847e33 +KPC_CM = 3.0856775814913673e21 +GYR_S = 3.15576e16 +G_KPC_KM2_S2_MSUN = 4.30091727003628e-6 +G_CGS = 6.6743e-8 + +# x = r / r_s at which the NFW circular velocity peaks. +NFW_X_AT_V_MAX = 2.1625815870646098 + +# The calibrated constant C of the Yang et al. (2024) collapse-time formula. +COLLAPSE_CALIBRATION_C = 0.75 + + +def _nu_eff_km_s_from(rho_s, scale_radius_kpc): + """ + The effective velocity dispersion nu_eff = 0.64 * V_max of the halo's CDM + baseline NFW profile (Yang et al. 2024; Outmezguine et al. 2023), in km/s. + + ``rho_s`` is in Msun / kpc^3 and ``scale_radius_kpc`` in kpc. + """ + x = NFW_X_AT_V_MAX + mass_factor = np.log1p(x) - x / (1.0 + x) + v_max_km_s = np.sqrt( + 4.0 + * np.pi + * G_KPC_KM2_S2_MSUN + * rho_s + * scale_radius_kpc**2 + * mass_factor + / x + ) + return 0.64 * v_max_km_s + + +def _sigma_eff_over_m_from(sigma_over_m, velocity_exponent, velocity_ref, nu_eff_km_s): + """ + The effective cross section per unit mass / of + Outmezguine et al. (2023), for a power-law velocity dependence + + sigma(v) = sigma_over_m * (v / velocity_ref)^(-velocity_exponent) + + averaged over a Maxwell-Boltzmann relative-velocity distribution + f(v) ~ v^2 exp(-v^2 / (4 nu_eff^2)). The average is closed form: + + sigma_eff/m = (sigma/m) (velocity_ref / (2 nu_eff))^a Gamma(4 - a/2) / 6 + + and requires ``velocity_exponent < 8`` for the average to converge. At + ``velocity_exponent = 0`` this reduces to the constant cross section. + """ + a = float(velocity_exponent) + + if a == 0.0: + return sigma_over_m + + if a >= 8.0: + raise ValueError( + "The Maxwell-Boltzmann average diverges for a " + f"velocity_exponent >= 8 (input value {a})." + ) + + return ( + sigma_over_m + * (velocity_ref / (2.0 * nu_eff_km_s)) ** a + * math.gamma(4.0 - a / 2.0) + / 6.0 + ) + + +def _tau_from(sigma_eff_over_m, t_age, rho_s, scale_radius_kpc): + """ + The dimensionless gravothermal time tau = t_age / t_c, via the Yang et al. + (2024) collapse timescale of the CDM baseline halo: + + t_c = (150 / C) / [ (sigma_eff/m) rho_s r_s sqrt(4 pi G rho_s) ] + + with C = 0.75. ``sigma_eff_over_m`` is in cm^2/g, ``t_age`` in Gyr, + ``rho_s`` in Msun / kpc^3 and ``scale_radius_kpc`` in kpc. + """ + if sigma_eff_over_m <= 0.0 or t_age <= 0.0: + return 0.0 + + rho_s_g_cm3 = rho_s * MSUN_G / KPC_CM**3 + scale_radius_cm = scale_radius_kpc * KPC_CM + + collapse_rate = ( + sigma_eff_over_m + * rho_s_g_cm3 + * scale_radius_cm + * np.sqrt(4.0 * np.pi * G_CGS * rho_s_g_cm3) + ) + t_c_seconds = (150.0 / COLLAPSE_CALIBRATION_C) / collapse_rate + + return float(t_age * GYR_S / t_c_seconds) + + +class YangSIDMMCRLudlowSph(YangSIDMSph): + r""" + The Yang et al. (2024) parametric SIDM halo constructed from physical + inputs, with the CDM baseline NFW parameters set by the Ludlow et al. + (2016) mass-concentration relation. + + The dimensionless gravothermal time is solved per halo as + :math:`\tau = t_{\rm age} / t_c` with + + .. math:: + t_c = \frac{150 / C}{(\sigma_{\rm eff}/m) \rho_s r_s + \sqrt{4 \pi G \rho_s}}, \quad C = 0.75, + + where :math:`(\rho_s, r_s)` are the halo's CDM baseline NFW parameters in + physical units. The effective cross section + :math:`\sigma_{\rm eff}/m` follows Outmezguine et al. (2023): a power-law + velocity dependence :math:`\sigma \propto v^{-a}` averaged over a + Maxwell-Boltzmann distribution at :math:`\nu_{\rm eff} = 0.64 V_{\rm max}`. + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + mass_at_200 + The mass of the halo enclosing 200 times the cosmic average density, + in solar masses. + sigma_over_m + The self-interaction cross section per unit mass, in cm^2 / g, + evaluated at the reference velocity ``velocity_ref``. + velocity_exponent + The power-law index a of the velocity-dependent cross section + sigma(v) = sigma_over_m * (v / velocity_ref)^(-a). Zero gives a + velocity-independent cross section (must be < 8 for the + Maxwell-Boltzmann average to converge). + velocity_ref + The reference velocity, in km/s, at which ``sigma_over_m`` is quoted. + Irrelevant when ``velocity_exponent`` is zero. + t_age + The age of the halo in Gyr, from which tau = t_age / t_c is computed. + For multi-plane systems construct each halo with the age appropriate to + its redshift plane. + redshift_object + The halo redshift, setting the Ludlow et al. (2016) concentration and + the arc-second unit conversions. + redshift_source + The source redshift, setting the critical surface density. + """ + + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + mass_at_200: float = 1e9, + sigma_over_m: float = 1.0, + velocity_exponent: float = 0.0, + velocity_ref: float = 10.0, + t_age: float = 10.0, + redshift_object: float = 0.5, + redshift_source: float = 1.0, + ): + self.mass_at_200 = mass_at_200 + self.sigma_over_m = sigma_over_m + self.velocity_exponent = velocity_exponent + self.velocity_ref = velocity_ref + self.t_age = t_age + self.redshift_object = redshift_object + self.redshift_source = redshift_source + + ( + concentration, + cosmic_average_density, + critical_surface_density, + kpc_per_arcsec, + ) = mcr_util.ludlow16_cosmology( + mass_at_200=mass_at_200, + redshift_object=redshift_object, + redshift_source=redshift_source, + xp=np, + ) + + radius_at_200 = ( + mass_at_200 / (200.0 * cosmic_average_density * (4.0 * np.pi / 3.0)) + ) ** (1.0 / 3.0) + scale_radius_kpc = radius_at_200 / concentration + + de_c = ( + 200.0 + / 3.0 + * ( + concentration**3 + / (np.log(1.0 + concentration) - concentration / (1.0 + concentration)) + ) + ) + rho_s = cosmic_average_density * de_c + kappa_s = rho_s * scale_radius_kpc / critical_surface_density + scale_radius = scale_radius_kpc / kpc_per_arcsec + + self.nu_eff = float( + _nu_eff_km_s_from(rho_s=rho_s, scale_radius_kpc=scale_radius_kpc) + ) + self.sigma_eff_over_m = float( + _sigma_eff_over_m_from( + sigma_over_m=sigma_over_m, + velocity_exponent=velocity_exponent, + velocity_ref=velocity_ref, + nu_eff_km_s=self.nu_eff, + ) + ) + + # The un-truncated tau (the base class truncates into the fitted + # [0, 1] range), kept for diagnostics such as core-collapse fractions. + self.tau_physical = _tau_from( + sigma_eff_over_m=self.sigma_eff_over_m, + t_age=t_age, + rho_s=rho_s, + scale_radius_kpc=scale_radius_kpc, + ) + tau = self.tau_physical + + super().__init__( + centre=centre, + kappa_s=kappa_s, + scale_radius=scale_radius, + tau=tau, + ) diff --git a/test_autogalaxy/profiles/mass/dark/test_yang24.py b/test_autogalaxy/profiles/mass/dark/test_yang24.py new file mode 100644 index 00000000..d4d5ea80 --- /dev/null +++ b/test_autogalaxy/profiles/mass/dark/test_yang24.py @@ -0,0 +1,257 @@ +import importlib.util + +import numpy as np +import pytest + +import autogalaxy as ag + +from autogalaxy.profiles.mass.dark.yang24 import _yang24_parameter_ratios_from +from autogalaxy.profiles.mass.dark.yang24_mcr import _sigma_eff_over_m_from + +# The `vmapped_deflections_from` path is jax-backed; these tests need jax +# installed to run (it ships via the `[optional]` extras). The NumPy-only +# Python-version matrix has no jax, so skip there rather than fail. +requires_jax = pytest.mark.skipif( + importlib.util.find_spec("jax") is None, + reason="requires jax (installed via the [optional] extras; absent on the NumPy-only matrix env)", +) + + +def test__parameter_ratios_are_exact_nfw_at_zero_tau(): + rho_s_ratio, r_s_ratio, r_c_ratio = _yang24_parameter_ratios_from(tau=0.0) + + assert rho_s_ratio == pytest.approx(1.0, abs=1.0e-12) + assert r_s_ratio == pytest.approx(1.0, abs=1.0e-12) + assert r_c_ratio == 0.0 + + +def test__zero_tau_limit_matches_nfw(): + yang = ag.mp.YangSIDMSph( + centre=(0.1, -0.2), + kappa_s=0.2, + scale_radius=2.0, + tau=0.0, + ) + nfw = ag.mp.NFWSph(centre=(0.1, -0.2), kappa_s=0.2, scale_radius=2.0) + + grid = ag.Grid2DIrregular([[0.5, 0.2], [1.0, -0.2], [2.0, 1.0]]) + + assert yang.convergence_2d_from(grid=grid) == pytest.approx( + nfw.convergence_2d_from(grid=grid).array, + rel=1.0e-8, + ) + assert yang.deflections_yx_2d_from(grid=grid) == pytest.approx( + nfw.deflections_yx_2d_from(grid=grid).array, + rel=1.0e-8, + ) + + +def test__evolved_halo_has_constant_density_core(): + profile = ag.mp.YangSIDMSph( + kappa_s=0.2, + scale_radius=2.0, + tau=0.5, + ) + + assert profile.core_radius_evolved > 0.0 + + density_centre = profile.density_3d_func(ag.ArrayIrregular([1.0e-6]))[0] + + # At the centre the beta = 4 cored form tends to the finite value + # rho_s * r_s / r_c. + analytic_centre = ( + profile.rho_s_evolved + * profile.scale_radius_evolved + / profile.core_radius_evolved + ) + + assert np.isfinite(density_centre) + assert density_centre == pytest.approx(analytic_centre, rel=1.0e-4) + + nfw_density = ag.mp.NFWSph(kappa_s=0.2, scale_radius=2.0).density_3d_func( + ag.ArrayIrregular([1.0e-6]) + )[0] + assert density_centre < nfw_density + + +def test__tau_is_truncated_into_fitted_range(): + profile = ag.mp.YangSIDMSph(kappa_s=0.2, scale_radius=2.0, tau=1.7) + truncated = ag.mp.YangSIDMSph(kappa_s=0.2, scale_radius=2.0, tau=1.0) + + assert profile.tau == 1.0 + assert profile.rho_s_evolved == truncated.rho_s_evolved + + assert ag.mp.YangSIDMSph(kappa_s=0.2, scale_radius=2.0, tau=-0.3).tau == 0.0 + + +def test__lensing_quantities_are_finite_and_positive(): + profile = ag.mp.YangSIDMSph( + kappa_s=0.2, + scale_radius=2.0, + tau=0.5, + ) + grid = ag.Grid2DIrregular([[0.0, 0.0], [0.2, 0.0], [1.0, 0.0]]) + + convergence = np.asarray(profile.convergence_2d_from(grid=grid).array) + deflections = np.asarray(profile.deflections_yx_2d_from(grid=grid).array) + potential = np.asarray(profile.potential_2d_from(grid=grid).array) + + assert np.isfinite(convergence).all() + assert np.isfinite(deflections).all() + assert np.isfinite(potential).all() + assert (convergence > 0.0).all() + assert deflections[0] == pytest.approx(np.array([0.0, 0.0]), abs=1.0e-8) + + +@requires_jax +def test__vmapped_deflections_match_instance_path_for_zero_tau(): + profile = ag.mp.YangSIDMSph( + centre=(0.1, -0.2), + kappa_s=0.2, + scale_radius=2.0, + tau=0.0, + ) + grid = np.array([[0.5, 0.2], [1.0, -0.2], [2.0, 1.0]]) + + params = np.array( + [ + [ + profile.centre[0], + profile.centre[1], + profile.rho_s_evolved, + profile.scale_radius_evolved, + profile.core_radius_evolved, + ] + ] + ) + mask = np.array([True]) + + vmapped = ag.mp.YangSIDMSph.vmapped_deflections_from( + grid=grid, + params_batch=params, + mask=mask, + ) + + np.testing.assert_allclose( + np.asarray(vmapped), + profile.deflections_yx_2d_from(grid=ag.Grid2DIrregular(grid)).array, + rtol=1.0e-6, + atol=1.0e-8, + ) + + +@requires_jax +def test__vmapped_deflections_match_instance_path_for_evolved_halo(): + profile = ag.mp.YangSIDMSph( + centre=(0.1, -0.2), + kappa_s=0.2, + scale_radius=2.0, + tau=0.5, + ) + grid = np.array([[0.5, 0.2], [1.0, -0.2], [2.0, 1.0]]) + + params = np.array( + [ + [ + profile.centre[0], + profile.centre[1], + profile.rho_s_evolved, + profile.scale_radius_evolved, + profile.core_radius_evolved, + ] + ] + ) + mask = np.array([True]) + + vmapped = ag.mp.YangSIDMSph.vmapped_deflections_from( + grid=grid, + params_batch=params, + mask=mask, + ) + + np.testing.assert_allclose( + np.asarray(vmapped), + profile.deflections_yx_2d_from(grid=ag.Grid2DIrregular(grid)).array, + rtol=5.0e-2, + atol=1.0e-3, + ) + + +def test__mcr_constructor_reduces_to_nfw_when_cross_section_is_zero(): + yang = ag.mp.YangSIDMMCRLudlowSph( + centre=(1.0, 2.0), + mass_at_200=1.0e9, + sigma_over_m=0.0, + t_age=10.0, + redshift_object=0.6, + redshift_source=2.5, + ) + nfw = ag.mp.NFWSph( + centre=(1.0, 2.0), + kappa_s=yang.kappa_s, + scale_radius=yang.scale_radius, + ) + + grid = ag.Grid2DIrregular([[1.0, 1.0], [2.0, 2.0]]) + + assert yang.tau == 0.0 + assert yang.deflections_yx_2d_from(grid=grid) == pytest.approx( + nfw.deflections_yx_2d_from(grid=grid).array, + rel=1.0e-8, + ) + + +def test__mcr_constructor_tau_increases_with_cross_section(): + taus = [ + ag.mp.YangSIDMMCRLudlowSph( + mass_at_200=1.0e9, + sigma_over_m=sigma_over_m, + t_age=10.0, + redshift_object=0.6, + redshift_source=2.5, + ).tau_physical + for sigma_over_m in [1.0, 10.0, 100.0] + ] + + assert taus[0] > 0.0 + assert taus[0] < taus[1] < taus[2] + assert taus[1] == pytest.approx(10.0 * taus[0], rel=1.0e-8) + + +def test__effective_cross_section_velocity_dependence(): + nu_eff = 20.0 + + constant = _sigma_eff_over_m_from( + sigma_over_m=3.0, + velocity_exponent=0.0, + velocity_ref=10.0, + nu_eff_km_s=nu_eff, + ) + assert constant == 3.0 + + # For sigma ~ v^-a quoted at a reference velocity below the halo's + # characteristic velocity, the effective cross section is suppressed. + suppressed = _sigma_eff_over_m_from( + sigma_over_m=3.0, + velocity_exponent=1.0, + velocity_ref=10.0, + nu_eff_km_s=nu_eff, + ) + assert suppressed < constant + + # An a = 2 power law has the closed form (v0 / 2 nu)^2 * Gamma(3) / 6. + analytic = 3.0 * (10.0 / (2.0 * nu_eff)) ** 2 * 2.0 / 6.0 + assert _sigma_eff_over_m_from( + sigma_over_m=3.0, + velocity_exponent=2.0, + velocity_ref=10.0, + nu_eff_km_s=nu_eff, + ) == pytest.approx(analytic, rel=1.0e-12) + + with pytest.raises(ValueError): + _sigma_eff_over_m_from( + sigma_over_m=3.0, + velocity_exponent=8.0, + velocity_ref=10.0, + nu_eff_km_s=nu_eff, + )