added kappa_s_and_scale_radius_and_core_radius

Author: Niek Wielders

autogalaxy/profiles/mass/dark/mcr_util.py

@@ -58,6 +58,56 @@ def kappa_s_and_scale_radius_for_duffy(mass_at_200, redshift_object, redshift_so
 
     return kappa_s, scale_radius, radius_at_200
 
+def kappa_s_and_scale_radius_and_core_radius(mass_at_200, f_c, redshift_object, redshift_source):
+    """
+    Computes the AutoGalaxy cNFW parameters (kappa_s, scale_radius, core_radius) for an 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
+    """
+
+    from astropy import units
+
+    from autogalaxy.cosmology.wrap import Planck15
+
+    cosmology = Planck15()
+
+    cosmic_average_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)
+
+    radius_at_200 = (
+        mass_at_200 / (200.0 * cosmic_average_density * (4.0 * np.pi / 3.0))
+    ) ** (
+        1.0 / 3.0
+    )  # r200
+    coefficient = 5.71 * (1.0 + redshift_object) ** (
+        -0.47
+    )  # The coefficient of Duffy mass-concentration (Duffy+2008)
+    concentration = coefficient * (mass_at_200 / 2.952465309e12) ** (
+        -0.084
+    )  # mass-concentration relation. (Duffy+2008)
+    mass_concentration_relation = ((f_c**2 * np.log(1 + concentration / f_c) + (1 - 2 * f_c) * np.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 * np.pi * scale_radius_kpc**3 * mass_concentration_relation)
+    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
+
 
 def _ludlow16_cosmology_callback(
     mass_at_200,