diff --git a/autogalaxy/config/priors/mass/dark/cnfw.yaml b/autogalaxy/config/priors/mass/dark/cnfw.yaml new file mode 100644 index 000000000..265c9782d --- /dev/null +++ b/autogalaxy/config/priors/mass/dark/cnfw.yaml @@ -0,0 +1,51 @@ +cNFWSph: + 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 + core_radius: + type: Uniform + lower_limit: 0.0 + upper_limit: 15.0 + width_modifier: + type: Relative + value: 0.2 + limits: + lower: 0.0 + upper: inf \ No newline at end of file diff --git a/autogalaxy/config/priors/mass/dark/cnfw_mcr.yaml b/autogalaxy/config/priors/mass/dark/cnfw_mcr.yaml new file mode 100644 index 000000000..4fd52d85d --- /dev/null +++ b/autogalaxy/config/priors/mass/dark/cnfw_mcr.yaml @@ -0,0 +1,61 @@ +cNFWMCRLudlowSph: + 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 + f_c: + type: Uniform + lower_limit: 0.0001 + upper_limit: 0.5 + width_modifier: + type: Relative + value: 0.2 + limits: + lower: 0.0001 + 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 \ No newline at end of file diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py index a4b358d0d..8ec0fa820 100644 --- a/autogalaxy/profiles/mass/__init__.py +++ b/autogalaxy/profiles/mass/__init__.py @@ -36,6 +36,8 @@ NFWMCRLudlow, gNFWMCRLudlow, NFWVirialMassConcSph, + cNFWSph, + cNFWMCRLudlowSph, ) from .stellar import ( Gaussian, diff --git a/autogalaxy/profiles/mass/dark/__init__.py b/autogalaxy/profiles/mass/dark/__init__.py index 37c26c344..5a7b29895 100644 --- a/autogalaxy/profiles/mass/dark/__init__.py +++ b/autogalaxy/profiles/mass/dark/__init__.py @@ -9,3 +9,5 @@ from .nfw_truncated_mcr import NFWTruncatedMCRLudlowSph, NFWTruncatedMCRDuffySph from .nfw_truncated_mcr_scatter import NFWTruncatedMCRScatterLudlowSph from .nfw_virial_mass_conc import NFWVirialMassConcSph +from .cnfw import cNFWSph +from .cnfw_mcr import cNFWMCRLudlowSph diff --git a/autogalaxy/profiles/mass/dark/cnfw.py b/autogalaxy/profiles/mass/dark/cnfw.py new file mode 100644 index 000000000..c98ad451b --- /dev/null +++ b/autogalaxy/profiles/mass/dark/cnfw.py @@ -0,0 +1,175 @@ +import numpy as np + +from typing import Tuple + +from autogalaxy.profiles.mass.dark.abstract import AbstractgNFW + +import autoarray as aa + +class cNFWSph(AbstractgNFW): + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + kappa_s: float = 0.05, + scale_radius: float = 1.0, + core_radius: float = 0.01, + ): + """ + Represents a spherical cored NFW density distribution + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + kappa_s + The overall normalization of the dark matter halo \| + (kappa_s = (rho_0 * D_d * scale_radius)/lensing_critical_density) + scale_radius + The cored NFW scale radius `theta_s`, as an angle on the sky in arcseconds. + core_radius + The cored NFW core radius `theta_c`, as an angle on the sky in arcseconds. + """ + + super().__init__( + centre=centre, + ell_comps=(0.0, 0.0)) + + self.kappa_s = kappa_s + self.scale_radius = scale_radius + self.core_radius = core_radius + + + @aa.grid_dec.to_vector_yx + @aa.grid_dec.transform + def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): + """ + Calculate the deflection angles on a grid of (y,x) arc-second coordinates. + + The input grid of (y,x) coordinates are transformed to a coordinate system centred on the profile centre with + and rotated based on the position angle defined from its `ell_comps` (this is described fully below). + + The numerical backend can be selected via the ``xp`` argument, allowing this + method to be used with both NumPy and JAX (e.g. inside ``jax.jit``-compiled + code). This is described fully later in this example. + + Parameters + ---------- + grid + The grid of (y,x) arc-second coordinates the deflection angles are computed on. + xp + The numerical backend to use, either `numpy` or `jax.numpy`. + """ + theta = self.radial_grid_from(grid=grid, xp=xp, **kwargs).array + theta = xp.maximum(theta, 1e-8) + + factor = ( + 4.0 + * self.kappa_s + * self.scale_radius**2 + ) + + deflection_r = ( + factor + * (self.F_func(theta, self.scale_radius, xp=xp) - self.F_func(theta, self.core_radius, xp=xp) + - (self.scale_radius - self.core_radius) * self.dev_F_func(theta, self.scale_radius, xp=xp) + ) + / (theta * (self.scale_radius - self.core_radius)**2) + ) + + + return self._cartesian_grid_via_radial_from( + grid=grid, + radius=deflection_r, + xp=xp, + **kwargs, + ) + + def F_func(self, theta, radius, xp=np): + + F = theta * 0.0 + + # theta < radius + mask1 = (theta > 0) & (theta < radius) + + # theta > radius + mask2 = theta > radius + + F = xp.where( + mask1, + ( + radius / 2 * xp.log(2 * radius / theta) + - xp.sqrt(radius ** 2 - theta ** 2) + * xp.arctanh(xp.sqrt((radius - theta) / (radius + theta))) + ), + F, + ) + + F = xp.where( + mask2, + ( + radius / 2 * xp.log(2 * radius / theta) + + xp.sqrt(theta ** 2 - radius ** 2) + * xp.arctan(xp.sqrt((theta - radius) / (theta + radius))) + ), + F, + ) + + return 2 * radius * F + + def dev_F_func(self, theta, radius, xp=np): + + dev_F = theta * 0.0 + + mask1 = (theta > 0) & (theta < radius) + mask2 = theta == radius + mask3 = theta > radius + + dev_F = xp.where( + mask1, + ( + radius * xp.log(2 * radius / theta) + - (2 * radius ** 2 - theta ** 2) / xp.sqrt(radius ** 2 - theta ** 2) + * xp.arctanh(xp.sqrt((radius - theta) / (radius + theta))) + ), + dev_F, + ) + + dev_F = xp.where( + mask2, + radius * (xp.log(2) - 1 / 2), + dev_F, + ) + + dev_F = xp.where( + mask3, + ( + radius * xp.log(2 * radius / theta) + + (theta ** 2 - 2 * radius ** 2) / xp.sqrt(theta ** 2 - radius ** 2) + * xp.arctan(xp.sqrt((theta - radius) / (theta + radius))) + ), + dev_F, + ) + + return 2 * dev_F + + @aa.grid_dec.to_array + def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): + """ + Convergence (dimensionless surface mass density) for the cored NFW profile. + This is not yet implemented for `cNFWSph`. + """ + raise NotImplementedError( + "convergence_2d_from is not implemented for cNFWSph; a physical cNFW " + "convergence expression must be added before use." + ) + + @aa.grid_dec.to_array + def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): + """ + Lensing potential for the cored NFW profile. + This is not yet implemented for `cNFWSph`. + """ + raise NotImplementedError( + "potential_2d_from is not implemented for cNFWSph; a physical cNFW " + "potential expression must be added before use." + ) \ No newline at end of file diff --git a/autogalaxy/profiles/mass/dark/cnfw_mcr.py b/autogalaxy/profiles/mass/dark/cnfw_mcr.py new file mode 100644 index 000000000..1a22e1056 --- /dev/null +++ b/autogalaxy/profiles/mass/dark/cnfw_mcr.py @@ -0,0 +1,34 @@ +from typing import Tuple + +from autogalaxy.profiles.mass.dark.cnfw import cNFWSph + +from autogalaxy.profiles.mass.dark import mcr_util + +class cNFWMCRLudlowSph(cNFWSph): + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + mass_at_200: float = 1e9, + f_c = 0.01, + redshift_object: float = 0.5, + redshift_source: float = 1.0, + ): + self.mass_at_200 = mass_at_200 + self.f_c = f_c + self.redshift_object = redshift_object + self.redshift_source = redshift_source + + ( + kappa_s, + scale_radius, + core_radius, + radius_at_200, + ) = mcr_util.kappa_s_scale_radius_and_core_radius_for_ludlow( + mass_at_200=mass_at_200, + scatter_sigma=0.0, + f_c=f_c, + redshift_object=redshift_object, + redshift_source=redshift_source, + ) + + super().__init__(centre=centre, kappa_s=kappa_s, scale_radius=scale_radius, core_radius=core_radius) \ No newline at end of file diff --git a/autogalaxy/profiles/mass/dark/mcr_util.py b/autogalaxy/profiles/mass/dark/mcr_util.py index 64468c039..9e76bca17 100644 --- a/autogalaxy/profiles/mass/dark/mcr_util.py +++ b/autogalaxy/profiles/mass/dark/mcr_util.py @@ -203,3 +203,68 @@ def kappa_s_and_scale_radius_for_ludlow( scale_radius = scale_radius_kpc / kpc_per_arcsec return kappa_s, scale_radius, radius_at_200 + +def kappa_s_scale_radius_and_core_radius_for_ludlow(mass_at_200, scatter_sigma, f_c, redshift_object, redshift_source): + """ + Computes the AutoGalaxy cNFW parameters (kappa_s, scale_radius, core_radius) for a cored NFW halo of the given + mass, enforcing the Penarrubia '12 mass-concentration relation. + + Interprets mass as *`M_{200c}`*, not `M_{200m}`. + + f_c = core_radius / scale radius + """ + + if isinstance(mass_at_200, (float, np.ndarray, np.float64)): + xp = np + else: + xp = jnp + + # ------------------------------------ + # Cosmology + concentration (callback) + # ------------------------------------ + + if xp is np: + ( + concentration, + cosmic_average_density, + critical_surface_density, + kpc_per_arcsec, + ) = _ludlow16_cosmology_callback( + mass_at_200, + redshift_object, + redshift_source, + ) + else: + ( + concentration, + cosmic_average_density, + critical_surface_density, + kpc_per_arcsec, + ) = ludlow16_cosmology_jax( + mass_at_200, + redshift_object, + redshift_source, + ) + + # Apply scatter (JAX-safe) + concentration = 10.0 ** (xp.log10(concentration) + scatter_sigma * 0.15) + + # ------------------------------------ + # JAX-native algebra + # ------------------------------------ + radius_at_200 = ( + mass_at_200 / (200.0 * cosmic_average_density * (4.0 * xp.pi / 3.0)) + ) ** ( + 1.0 / 3.0 + ) # r200 + + mcr_penarrubia = ((f_c**2 * xp.log(1 + concentration / f_c) + (1 - 2 * f_c) * xp.log(1 + concentration)) / (1 + f_c)**2 + - concentration / ((1+concentration) * (1-f_c))) #mass concentration relation (Penarrubia+2012) + + scale_radius_kpc = radius_at_200 / concentration # scale radius in kpc + rho_0 = mass_at_200 / (4 * xp.pi * scale_radius_kpc**3 * mcr_penarrubia) + kappa_s = rho_0 * scale_radius_kpc / critical_surface_density # kappa_s + scale_radius = scale_radius_kpc / kpc_per_arcsec # scale radius in arcsec + core_radius = f_c * scale_radius # core radius in arcsec + + return kappa_s, scale_radius, core_radius, radius_at_200 diff --git a/test_autogalaxy/profiles/mass/dark/test_cnfw.py b/test_autogalaxy/profiles/mass/dark/test_cnfw.py new file mode 100644 index 000000000..2dae506c8 --- /dev/null +++ b/test_autogalaxy/profiles/mass/dark/test_cnfw.py @@ -0,0 +1,12 @@ +import numpy as np +import pytest + +import autogalaxy as ag + +def test__deflections_yx_2d_from(): + cnfw = ag.mp.cNFWSph(centre=(0.0, 0.0), kappa_s=0.01591814312464436, scale_radius=0.36, core_radius=0.036) + + deflection_2d = cnfw.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[1.0, 0.0]])) + deflection_r = np.sqrt(deflection_2d[0, 0]**2 + deflection_2d[0, 1]**2) + + assert deflection_r == pytest.approx(0.006034319441107217, 1.0e-8) diff --git a/test_autogalaxy/profiles/mass/dark/test_nfw_mcr.py b/test_autogalaxy/profiles/mass/dark/test_nfw_mcr.py index 143f6f108..47964eb69 100644 --- a/test_autogalaxy/profiles/mass/dark/test_nfw_mcr.py +++ b/test_autogalaxy/profiles/mass/dark/test_nfw_mcr.py @@ -38,7 +38,7 @@ def test__mass_and_concentration_consistent_with_normal_nfw(): redshift_profile=0.6, redshift_source=2.5, cosmology=cosmology ) - # We uare using the NFWTruncatedSph to check the mass gives a conosistnt kappa_s, given certain radii. + # We are using the NFWTruncatedSph to check that the mass gives a consistent kappa_s, given certain radii. assert mass_at_200_via_kappa_s == mass_at_200_via_mass assert concentration_via_kappa_s == concentration_via_mass @@ -216,3 +216,58 @@ def test__same_as_above_but_generalized_elliptical(): deflections = nfw_kappa_s.deflections_yx_2d_from(grid=grid) assert (deflections_ludlow == deflections).all() + +def test__same_as_above_but_cored_nfw(): + + from autogalaxy.cosmology.model import FlatLambdaCDMWrap + + cosmology = FlatLambdaCDMWrap(H0=70.0, Om0=0.3) + + mp = ag.mp.cNFWMCRLudlowSph( + centre=(1.0, 2.0), + mass_at_200=1.0e9, + f_c=0.01, + redshift_object=0.6, + redshift_source=2.5, + ) + + mass_at_200_via_mass = mp.mass_at_200_solar_masses( + redshift_object=0.6, redshift_source=2.5, cosmology=cosmology + ) + concentration_via_mass = mp.concentration( + redshift_profile=0.6, redshift_source=2.5, cosmology=cosmology + ) + + cnfw_kappa_s = ag.mp.cNFWSph( + centre=(1.0, 2.0), + kappa_s=mp.kappa_s, + scale_radius=mp.scale_radius, + core_radius=mp.core_radius, + ) + + mass_at_200_via_kappa_s = cnfw_kappa_s.mass_at_200_solar_masses( + redshift_object=0.6, redshift_source=2.5, cosmology=cosmology + ) + concentration_via_kappa_s = cnfw_kappa_s.concentration( + redshift_profile=0.6, redshift_source=2.5, cosmology=cosmology + ) + + # We are using the NFWTruncatedSph to check the mass gives a consistent kappa_s, given certain radii. + + assert mass_at_200_via_kappa_s == mass_at_200_via_mass + assert concentration_via_kappa_s == concentration_via_mass + + assert mp.centre == (1.0, 2.0) + + assert mp.axis_ratio() == 1.0 + + assert mp.angle() == 0.0 + + assert mp.scale_radius == pytest.approx(0.21158, 1.0e-4) + + assert mp.core_radius == pytest.approx(0.0021158, 1.0e-4) + + deflections_ludlow = mp.deflections_yx_2d_from(grid=grid) + deflections = cnfw_kappa_s.deflections_yx_2d_from(grid=grid) + + assert (deflections_ludlow == deflections).all()