Merge pull request #190 from Jammy2211/feature/kassidy_nfw

Feature/kassidy nfw

autogalaxy/config/notation.yaml

@@ -8,6 +8,7 @@ label:
     centre_0: y
     centre_1: x
     coefficient: \lambda
+    c_2: c_{\rm 2}
     concentration: conc
     core_radius: C_{\rm r}
     core_radius_0: C_{rm r0}
@@ -27,6 +28,7 @@ label:
     intensity: I_{\rm b}
     kappa: \kappa
     kappa_s: \kappa_{\rm s}
+    log10m_vir: log_{\rm 10}(m_{vir})
     m: m
     mass: M
     mass_at_200: M_{\rm 200}
@@ -41,6 +43,7 @@ label:
     normalization_scale: n
     outer_coefficient: \lambda_{\rm 2}
     outer_slope: t_{\rm 2}
+    overdens: \Delta_{\rm vir}
     pixels: N_{\rm pix}
     radius_break: R_{\rm b}
     redshift: z
@@ -77,6 +80,7 @@ label_format:
     angle_binary: '{:.4f}'
     angular_diameter_distance_to_earth: '{:.4f}'
     beta: '{:.4f}'
+    c_2: '{:.4f}'
     centre_0: '{:.4f}'
     centre_1: '{:.4f}'
     coefficient: '{:.4f}'
@@ -99,6 +103,7 @@ label_format:
     kappa: '{:.4f}'
     kappa_s: '{:.4f}'
     kpc_per_arcsec: '{:.4f}'
+    log10m_vir: '{:.4f}'
     luminosity: '{:.4e}'
     m: '{:.1f}'
     mass: '{:.4e}'
@@ -115,6 +120,7 @@ label_format:
     normalization_scale: '{:.4f}'
     outer_coefficient: '{:.4f}'
     outer_slope: '{:.4f}'
+    overdens: '{:.4f}'
     pixels: '{:.4f}'
     radius: '{:.4f}'
     radius_break: '{:.4f}'

autogalaxy/config/priors/mass/dark/gnfw_virial_mass_conc.yaml

@@ -0,0 +1,85 @@
+gNFWVirialMassConcSph:
+  centre_0:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    lower_limit: -inf
+    upper_limit: inf
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    gaussian_limits:
+      lower: -inf
+      upper: inf
+  centre_1:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    lower_limit: -inf
+    upper_limit: inf
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    gaussian_limits:
+      lower: -inf
+      upper: inf
+  log10m_vir:
+    type: Uniform
+    lower_limit: 7.0
+    upper_limit: 12.0
+    width_modifier:
+      type: Relative
+      value: 0.5
+    gaussian_limits:
+      lower: 0.0
+      upper: inf
+  c_2:
+    type: LogUniform
+    lower_limit: 1.0
+    upper_limit: 100.0
+    width_modifier:
+      type: Relative
+      value: 0.5
+    gaussian_limits:
+      lower: 0.0
+      upper: inf
+  overdens:
+    type: Uniform
+    lower_limit: 100.0
+    upper_limit: 250.0
+    width_modifier:
+      type: Relative
+      value: 0.5
+    gaussian_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
+    gaussian_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
+    gaussian_limits:
+      lower: 0.0
+      upper: inf
+  inner_slope:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 2.0
+    width_modifier:
+      type: Absolute
+      value: 0.3
+    gaussian_limits:
+      lower: -1.0
+      upper: 3.0

autogalaxy/profiles/light/snr/__init__.py

@@ -11,4 +11,4 @@
     ElsonFreeFall,
     ElsonFreeFallSph,
 )
-from .sersic_core import SersicCore
+from .sersic_core import SersicCore

autogalaxy/profiles/mass/__init__.py

@@ -16,6 +16,7 @@
 from .dark import (
     gNFW,
     gNFWSph,
+    gNFWVirialMassConcSph,
     NFWTruncatedSph,
     NFWTruncatedMCRDuffySph,
     NFWTruncatedMCRLudlowSph,

autogalaxy/profiles/mass/dark/__init__.py

@@ -1,5 +1,6 @@
 from .gnfw import gNFW, gNFWSph
 from .gnfw_mcr import gNFWMCRLudlow
+from .gnfw_virial_mass_conc import gNFWVirialMassConcSph
 from .nfw import NFW, NFWSph
 from .nfw_mcr import NFWMCRLudlowSph, NFWMCRDuffySph, NFWMCRLudlow
 from .nfw_mcr_scatter import NFWMCRScatterLudlow, NFWMCRScatterLudlowSph

autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py

@@ -0,0 +1,133 @@
+from typing import Tuple
+
+from autogalaxy.profiles.mass.dark.gnfw import gNFWSph
+
+from astropy import units
+
+import numpy as np
+from autogalaxy import cosmology as cosmo
+
+from scipy.integrate import quad
+
+
+def kappa_s_and_scale_radius(
+    cosmology, virial_mass, c_2, overdens, redshift_object, redshift_source, inner_slope
+):
+    concentration = (2 - inner_slope) * c_2  # gNFW concentration
+
+    critical_density = (
+        cosmology.critical_density(redshift_object).to(units.solMass / units.kpc**3)
+    ).value
+
+    critical_surface_density = (
+        cosmology.critical_surface_density_between_redshifts_solar_mass_per_kpc2_from(
+            redshift_0=redshift_object, redshift_1=redshift_source
+        )
+    )
+
+    kpc_per_arcsec = cosmology.kpc_per_arcsec_from(redshift=redshift_object)
+
+    if overdens == 0:
+        x = cosmology.Om(redshift_object) - 1
+        overdens = 18 * np.pi**2 + 82 * x - 39 * x**2  # Bryan & Norman (1998)
+
+    virial_radius = (
+        virial_mass / (overdens * critical_density * (4.0 * np.pi / 3.0))
+    ) ** (
+        1.0 / 3.0
+    )  # r_vir
+
+    scale_radius_kpc = (
+        virial_radius / concentration
+    )  # scale radius of gNFW profile in kpc
+
+    ##############################
+    def integrand(r):
+        return (r**2 / r**inner_slope) * (1 + r / scale_radius_kpc) ** (
+            inner_slope - 3
+        )
+
+    de_c = (
+        (overdens / 3.0)
+        * (virial_radius**3 / scale_radius_kpc**inner_slope)
+        / quad(integrand, 0, virial_radius)[0]
+    )  # rho_c
+    ##############################
+
+    rho_s = critical_density * de_c  # rho_s
+    kappa_s = rho_s * scale_radius_kpc / critical_surface_density  # kappa_s
+    scale_radius = scale_radius_kpc / kpc_per_arcsec  # scale radius in arcsec
+
+    return kappa_s, scale_radius, virial_radius, overdens
+
+
+class gNFWVirialMassConcSph(gNFWSph):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        log10m_vir: float = 12.0,
+        c_2: float = 10.0,
+        overdens: float = 0.0,
+        redshift_object: float = 0.5,
+        redshift_source: float = 1.0,
+        inner_slope: float = 1.0,
+    ):
+        """
+        Spherical gNFW profile initialized with the virial mass and c_2 concentration of the halo.
+
+        The virial radius of the halo is defined as the radius at which the density of the halo
+        equals overdens * the critical density of the Universe. r_vir = (3*m_vir/4*pi*overdens*critical_density)^1/3.
+
+        If the overdens parameter is set to 0, the virial overdensity of Bryan & Norman (1998) will be used.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) arc-second coordinates of the profile centre.
+        log10m_vir
+            The log10(virial mass) of the dark matter halo.
+        c_2
+            The c_2 concentration of the dark matter halo, which equals r_vir/r_2, where r_2 is the
+            radius at which the logarithmic density slope equals -2.
+        overdens
+            The spherical overdensity used to define the virial radius of the dark matter
+            halo: r_vir = (3*m_vir/4*pi*overdens*critical_density)^1/3. If this parameter is set to 0, the virial
+            overdensity of Bryan & Norman (1998) will be used.
+        redshift_object
+            Lens redshift.
+        redshift_source
+            Source redshift.
+        inner_slope
+            The inner slope of the dark matter halo's gNFW density profile.
+        """
+
+        self.log10m_vir = log10m_vir
+        self.c_2 = c_2
+        self.redshift_object = redshift_object
+        self.redshift_source = redshift_source
+        self.inner_slope = inner_slope
+
+        (
+            kappa_s,
+            scale_radius,
+            virial_radius,
+            overdens,
+        ) = kappa_s_and_scale_radius(
+            cosmology=cosmo.Planck15(),
+            virial_mass=10**log10m_vir,
+            c_2=c_2,
+            overdens=overdens,
+            redshift_object=redshift_object,
+            redshift_source=redshift_source,
+            inner_slope=inner_slope,
+        )
+
+        self.virial_radius = virial_radius
+        self.overdens = overdens
+
+        super().__init__(
+            centre=centre,
+            kappa_s=kappa_s,
+            inner_slope=inner_slope,
+            scale_radius=scale_radius,
+        )

test_autogalaxy/config/notation.yaml

@@ -25,6 +25,7 @@ label:
     intensity: I_{\rm b}
     kappa: \kappa
     kappa_s: \kappa_{\rm s}
+    log10m_vir: log_{\rm 10}(m_{vir})
     m: m
     mass: M
     mass_at_200: M_{\rm 200}
@@ -39,6 +40,7 @@ label:
     normalization_scale: n
     outer_coefficient: \lambda_{\rm 2}
     outer_slope: t_{\rm 2}
+    overdens: \Delta_{\rm vir}
     pixels: N_{\rm pix}
     radius_break: R_{\rm b}
     redshift: z

test_autogalaxy/profiles/mass/dark/test_gnfw_virial_mass_conc.py

@@ -0,0 +1,25 @@
+import numpy as np
+import pytest
+
+import autogalaxy as ag
+
+grid = ag.Grid2DIrregular([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [2.0, 4.0]])
+
+
+def test__deflections_2d_via_integral_from():
+    mp = ag.mp.gNFWVirialMassConcSph(
+        centre=(0.0, 0.0),
+        log10m_vir=12.0,
+        c_2=10.0,
+        overdens=0.0,
+        redshift_object=0.5,
+        redshift_source=1.0,
+        inner_slope=1.0,
+    )
+
+    deflections = mp.deflections_2d_via_integral_from(
+        grid=ag.Grid2DIrregular([[0.1875, 0.1625]])
+    )
+
+    assert deflections[0, 0] == pytest.approx(0.0466231, 1e-3)
+    assert deflections[0, 1] == pytest.approx(0.04040671, 1e-3)