diff --git a/clmm/model_profile.py b/clmm/model_profile.py
new file mode 100644
index 000000000..7ce24a1dc
--- /dev/null
+++ b/clmm/model_profile.py
@@ -0,0 +1,300 @@
+import numpy as np
+from astropy import units as un
+import clmm
+from clmm.utils.constants import Constants as const
+from clmm.dataops import _compute_tangential_shear, _compute_cross_shear
+
+
+class profiles:
+ '''
+ Class for computing shear profiles
+
+ Attributes
+ ----------
+
+ cosmo : clmm.cosmology.Cosmology object
+ CLMM Cosmology object
+ mdelta : float
+ Galaxy cluster mass in :math:`M_\odot`.
+ cdelta : float
+ Galaxy cluster concentration
+ z_l: float
+ Redshift of the lens
+ z_s: float
+ Redshift of the source
+ halo_profile_model : str, optional
+ Profile model parameterization (letter case independent):
+ * 'nfw' (default)
+ * 'einasto' - not in cluster_toolkit
+ * 'hernquist' - not in cluster_toolkit
+ delta_mdef : int, optional
+ Mass overdensity definition; defaults to 200.
+ massdef : str, optional
+ Profile mass definition, with the following supported options (letter case independent):
+ * 'mean' (default)
+ * 'critical'
+ * 'virial'
+ compute_sigma: bool
+ True for computing Sigma (surface density) profiles,
+ False for computing kappa (convergence) profiles
+ '''
+
+ def __init__(self,cosmo,mdelta,cdelta,z_l,z_s,halo_profile_model='nfw',delta_mdef=200,massdef='mean',compute_sigma=False):
+
+ self.mdelta = mdelta # M_sun
+ self.cdelta = cdelta
+ self.z_l = z_l
+ self.z_s = z_s
+ self.halo_profile_model = halo_profile_model
+ self.delta_mdef = delta_mdef
+ self.massdef = massdef
+ self.compute_sigma = compute_sigma
+
+ # physical constants with units (cosmological)
+ self.const_G = const.GNEWT_SOLAR_MASS* (un.Mpc/const.PC_TO_METER/1e6)**3 /un.M_sun / un.s**2 # Mpc^3/Msun/s^2
+ self.const_c = const.CLIGHT * (un.Mpc/const.PC_TO_METER/1e6)/ un.s # Mpc/s
+
+ # cosmology
+ self.cosmo = cosmo
+ return
+
+
+ def get_r_vals(self,nbins,lower_r,upper_r,log=False):
+ '''
+ Create radial bins for profiles
+
+ Attributes
+ ----------
+ nbins: int
+ number of bins
+ lower_r: float
+ lowest bin value
+ upper_r: float
+ highest bin value
+
+ Returns
+ ----------
+ r_vals: float or array
+ bin edges in units provided, or kpc if none are provided
+ '''
+
+ #create r values
+ self.nbins = nbins
+ self.lower_r = lower_r
+ self.upper_r = upper_r
+ self.log_r_vals = log
+
+ if self.log_r_vals:
+ if self.lower_r == 0:
+ raise ValueError('log(lower_r) cannot be evaluated for lower_r = 0')
+
+ r_vals = np.geomspace(self.lower_r,self.upper_r,self.nbins+1)
+
+ else:
+ r_vals = np.linspace(self.lower_r,self.upper_r,self.nbins+1)
+
+ #if no units are provided, set units to kpc, otherwise keep input units
+ r_vals *= un.dimensionless_unscaled
+ if r_vals.unit == un.dimensionless_unscaled:
+ r_vals *= un.kpc
+
+ self.r_vals = r_vals
+ return self.r_vals
+
+ def kappa_or_profile(self,r_vals=None):
+ '''
+ Compute kappa or surface density profiles
+
+ Attributes
+ ----------
+ r_vals: float or array
+ radial bin values for halo lensing profile
+
+ Returns
+ ----------
+ nfw: float or array
+ If compute_sigma is True, returns the surface density (Sigma) profile,
+ else, returns kappa profiles
+ '''
+
+ if np.all(r_vals.value) == None:
+ r_vals = self.r_vals
+
+ if self.compute_sigma:
+ return clmm.compute_surface_density(r_vals.value,self.mdelta.value,self.cdelta,
+ self.z_l,self.cosmo,self.delta_mdef,
+ self.halo_profile_model,self.massdef) * un.M_sun/un.Mpc**2
+ else:
+ return clmm.compute_convergence(r_vals.value,self.mdelta.value,self.cdelta,
+ self.z_l,self.z_s,self.cosmo,self.delta_mdef,
+ self.halo_profile_model,self.massdef) * un.dimensionless_unscaled
+
+
+
+ def make_grid(self,npix,r_max):
+ '''
+ Make 2D grid centered on origin
+
+ Attributes
+ ----------
+ npix: int
+ resolution of grid
+ r_max: float
+ width of grid is 2*r_max
+
+ Returns
+ ---------
+ r_2D: array
+ 2D grid (side=npix)
+ '''
+
+ npix *= 1j #imaginary so that ogrid knows to interpret npix as array size
+ y,x = np.ogrid[-r_max:r_max:npix,-r_max:r_max:npix]
+ r_2D = np.hypot(y,x)
+
+ return r_2D
+
+ def kappa_or_profile_2D(self,npix,r_max):
+ '''
+ Create 2D NFW profile on a grid
+
+ Attributes
+ ----------
+ npix: int
+ resolution of grid
+ r_max: float
+ width of grid is 2*r_max
+
+ Returns
+ ---------
+ kappa_or_profile: array
+ NFW profile (Sigma or kappa) for each point of the 2D grid
+
+ '''
+
+ self.npix = npix
+ self.r_max = r_max
+ r_2D = self.make_grid(self.npix,self.r_max)
+
+ kappa = np.array([self.kappa_or_profile(r_2D[i]) for i in range(self.npix)])
+ return kappa
+
+
+
+def KaiserSquires(Sigma):
+ '''
+ Kaiser & Squires method for computing shear components
+
+ Attributes
+ ----------
+ Sigma: array
+ Sigma map
+
+ Returns
+ ---------
+ e1,e2: arrays
+ shear maps
+ '''
+
+ kappa_tilde = np.fft.fft2(Sigma)
+
+ k = np.fft.fftfreq(kappa_tilde.shape[0])
+
+ oper_1 = - 1./(k[:, None]**2 + k[None, :]**2) * (k[:, None]**2 - k[None, :]**2)
+ oper_2 = - 2./(k[:, None]**2 + k[None, :]**2) * k[:, None]*k[None, :]
+
+ oper_1[0, 0] = -1
+ oper_2[0, 0] = -1
+
+ e1_tilde = oper_1*kappa_tilde
+ e2_tilde = oper_2*kappa_tilde
+
+ e1 = np.fft.ifft2(e1_tilde).real
+ e2 = np.fft.ifft2(e2_tilde).real
+
+ return e1, e2
+
+def get_Tangential_and_Cross(e1, e2, center, dx=10./1000.):
+ '''
+ Compute tangential shear profiles
+
+ Attributes
+ ----------
+ e1,e2: arrays
+ shear maps
+ center: array
+ center of the maps
+ dx: float
+ physcial size of a pixel
+
+ Returns
+ ---------
+ et,ex: arrays
+ tangential and cross shear maps aroud the given center
+ radius, angle: arrays
+ radius and angle maps, used for computing the radial profile
+ '''
+
+ n = e1.shape[0]
+ xx = np.arange(-n/2, n/2)*dx
+ XX, YY = np.meshgrid(xx, xx)
+
+ center_1, center_2 = center
+ from_cent_1 = XX - center_1
+ from_cent_2 = YY - center_2
+
+ angle = -np.sign(from_cent_2) * np.arccos(from_cent_1/np.sqrt(from_cent_1**2 + from_cent_2**2))
+ radius = np.sqrt(from_cent_1**2 + from_cent_2**2)
+
+ angle[np.isnan(angle)] = 0
+
+ et = _compute_tangential_shear(e1, e2, angle)
+ ex = _compute_cross_shear(e1, e2, angle)
+
+ return et, ex, radius, angle
+
+def get_Radial(radius_map,angle_map,r_bins,phi_bins,kappa_map,et_map):
+ '''
+ Calculate the radial convergence profile from the kappa map
+ and the radial tangential shear profile from the e_tangential map
+
+ Attributes
+ ----------
+ radius_map: array
+ radius map
+ angle_map:
+ angle map
+ r_bins: array
+ radial bins
+ phi_bins: array
+ angular bins
+ kappa_map: array
+ convergence map
+ et_map: array
+ tangential shear map
+
+ Returns
+ ---------
+ kappa_radial: array
+ radial convergence profile
+ gammat_radial: array
+ radial tangential shear profile
+ '''
+
+ _N = np.zeros((len(r_bins)-1, len(phi_bins)-1))
+ _K = np.zeros((len(r_bins)-1, len(phi_bins)-1))
+ _GT = np.zeros((len(r_bins)-1, len(phi_bins)-1))
+
+ ii_r = np.digitize(radius_map, r_bins) - 1
+ ii_phi = np.digitize(angle_map, phi_bins) - 1
+
+ for i_r, i_phi, kk, gg in zip(ii_r.flatten(), ii_phi.flatten(), kappa_map.flatten(), et_map.flatten()):
+
+ _N[i_r, i_phi] += 1
+ _K[i_r, i_phi] += kk
+ _GT[i_r, i_phi] += gg
+
+ kappa_radial = _K/_N
+ gammat_radial = _GT/_N
+
+ return kappa_radial,gammat_radial
diff --git a/examples/demo_model_profile.ipynb b/examples/demo_model_profile.ipynb
new file mode 100644
index 000000000..c37bdae21
--- /dev/null
+++ b/examples/demo_model_profile.ipynb
@@ -0,0 +1,17577 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "45308159",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib inline\n",
+ "%config InlineBackend.figure_format = 'svg'\n",
+ "\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import sys\n",
+ "from astropy import units as un\n",
+ "import os\n",
+ "\n",
+ "os.environ[\n",
+ " \"CLMM_MODELING_BACKEND\"\n",
+ "] = \"ccl\" # here you may choose ccl, nc (NumCosmo) or ct (cluster_toolkit)\n",
+ "\n",
+ "import clmm\n",
+ "\n",
+ "from clmm import Cosmology\n",
+ "from clmm.support import mock_data as mock\n",
+ "from clmm.galaxycluster import GalaxyCluster\n",
+ "from clmm import dataops\n",
+ "\n",
+ "# profile modelling package\n",
+ "import clmm.model_profile as mod"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7ba07a0d",
+ "metadata": {},
+ "source": [
+ "Make sure we know which version we're using"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "bbd94d3a",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "'0.16.0'"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "clmm.__version__"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "66c330d7",
+ "metadata": {},
+ "source": [
+ "The modelling module can be used to create 2D lensing maps. This allows us to forward model systematics, starting at the map level. In this notebook, we use miscentering as an example systematic."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 249,
+ "id": "909edd57",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# specify parameters for NFW halo\n",
+ "c = 3.5 #concentration\n",
+ "z_l = 0.5 #lens redshift\n",
+ "z_s = 1.0 #source redshift\n",
+ "Mass = 1e14 * un.M_sun\n",
+ "\n",
+ "# define physical size of grid\n",
+ "r_min = 1e-3 * un.Mpc\n",
+ "r_max = 3 * un.Mpc\n",
+ "\n",
+ "# number of pixels in grid\n",
+ "npix = 512\n",
+ "\n",
+ "# define physcial size of a pixel\n",
+ "dx = r_max/npix\n",
+ "\n",
+ "# define the center of the map (in Mpc)\n",
+ "center_1 = 0\n",
+ "center_2 = 0\n",
+ "\n",
+ "# define radial and angular bins\n",
+ "N_r_bins = 30\n",
+ "N_phi_bins = 25\n",
+ "r_bins = np.insert(np.logspace(-2, np.log10(r_max.value), N_r_bins), 0, r_min.value) * r_max.unit\n",
+ "phi_bins = np.linspace(-np.pi, np.pi+0.01, num=N_phi_bins+1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4e244831-d6ea-4e45-a0cf-1eb80979cb6b",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "id": "de3d76d6",
+ "metadata": {},
+ "source": [
+ "Now we can create the convergence maps"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 250,
+ "id": "5b51bbc8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# set cosmology\n",
+ "cosmo = Cosmology(H0=70.0, Omega_dm0=0.25, Omega_b0=0.05, Omega_k0=0.0)\n",
+ "\n",
+ "# lensing_map has the option to return either Sigma map or the convergence map\n",
+ "lensing_map = mod.profiles(cosmo,mdelta=Mass,cdelta=c,z_l=z_l,z_s=z_s,\n",
+ " halo_profile_model='nfw',delta_mdef=200,\n",
+ " massdef='critical',compute_sigma=False)\n",
+ "\n",
+ "# make map\n",
+ "kappa_map = lensing_map.kappa_or_profile_2D(npix,r_max) "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 265,
+ "id": "4dff16b1-28de-453f-bdbd-e5972eb434f8",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "axlim = r_max.value\n",
+ "plt.imshow(kappa_map,extent=[-axlim,axlim,-axlim,axlim],origin='lower')\n",
+ "plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ "plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ "plt.colorbar().set_label('$\\\\kappa$',rotation=0,fontsize=14)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e073090a",
+ "metadata": {},
+ "source": [
+ "Repeat the procedure for a Sigma map"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 252,
+ "id": "2aaf9cb0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# lensing_map has the option to return either Sigma map or the convergence map.\n",
+ "lensing_map = mod.profiles(cosmo,mdelta=Mass,cdelta=c,z_l=z_l,z_s=z_s,\n",
+ " halo_profile_model='nfw',delta_mdef=200,\n",
+ " massdef='critical',compute_sigma=True)\n",
+ "\n",
+ "# make map\n",
+ "Sigma_map = lensing_map.kappa_or_profile_2D(npix,r_max)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 253,
+ "id": "72993065",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "axlim = r_max.value\n",
+ "plt.imshow(Sigma_map,extent=[-axlim,axlim,-axlim,axlim],origin='lower')\n",
+ "plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ "plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ "plt.colorbar().set_label('$\\\\Sigma [M_\\odot\\; Mpc^{-2}]$',rotation=90)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4b28b017",
+ "metadata": {},
+ "source": [
+ "We can transform from the convergence to the shear with the Kaiser Squires algorithm, returning the e1 and e2 shear maps."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 254,
+ "id": "eec34ea3",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:230: RuntimeWarning: divide by zero encountered in divide\n",
+ " ---------\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:230: RuntimeWarning: invalid value encountered in multiply\n",
+ " ---------\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:231: RuntimeWarning: divide by zero encountered in divide\n",
+ " et,ex: arrays\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:231: RuntimeWarning: invalid value encountered in multiply\n",
+ " et,ex: arrays\n"
+ ]
+ }
+ ],
+ "source": [
+ "e1_map, e2_map = mod.KaiserSquires(kappa_map)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "25f42de1",
+ "metadata": {},
+ "source": [
+ "Plot the e1 and e2 maps"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 255,
+ "id": "2ba6da86",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.imshow(e1_map)\n",
+ "plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ "plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ "plt.colorbar().set_label('$e_1$',rotation=90)\n",
+ "plt.show()\n",
+ "\n",
+ "plt.imshow(e2_map)\n",
+ "plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ "plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ "plt.colorbar().set_label('$e_2$',rotation=90)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "12c04b61",
+ "metadata": {},
+ "source": [
+ "from the convergence map we can compute the convergence profile, and from the e1 and e2 maps we can compute the tangential shear profile."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f6af6b7e",
+ "metadata": {},
+ "source": [
+ "The getTangential function uses the e1 and e2 maps as an input, along with a map center, and returns e_tangential and e_cross maps about a given center. \n",
+ "\n",
+ "The function also returns a radius map and an angle map, which is used when computing radial profile with getRadial"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 256,
+ "id": "6b8c1825",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:273: RuntimeWarning: invalid value encountered in divide\n",
+ " tangential shear map\n"
+ ]
+ }
+ ],
+ "source": [
+ "et_map, ex_map, radius_map, angle_map = mod.get_Tangential_and_Cross(e1_map, e2_map, [center_1, center_2], dx=dx.value)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 257,
+ "id": "e743f1b8",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.imshow(et_map)\n",
+ "plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ "plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ "plt.colorbar().set_label('$e_t$',rotation=90)\n",
+ "plt.show()\n",
+ "\n",
+ "plt.imshow(ex_map)\n",
+ "plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ "plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ "plt.colorbar().set_label('$e_x$',rotation=90)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cde1c95f",
+ "metadata": {},
+ "source": [
+ "With the radius_map and angle_map, we can calculate the radial convergence profile from the kappa map and the radial tangential shear profile from the e_tangential map, using the getRadial function"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 258,
+ "id": "e3d066a0",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n"
+ ]
+ }
+ ],
+ "source": [
+ "kappa_radial, gammat_radial = mod.get_Radial(radius_map,angle_map,\n",
+ " r_bins.value,phi_bins,\n",
+ " kappa_map,et_map)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f89da817",
+ "metadata": {},
+ "source": [
+ "plot the radial convergence and tangential shear profiles. Note that these profiles are not averaged over all angles yet, meaning that these profiles are 2D (radial and angular). This is useful for miscentering which leaves anisotropic signatures on the resulting lensing profiles"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 259,
+ "id": "080ff2f4",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Text(0.5, 1.0, '$\\\\gamma_t$')"
+ ]
+ },
+ "execution_count": 259,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig,ax = plt.subplots(1,2,figsize=(9,4))\n",
+ "\n",
+ "ax[0].imshow(kappa_radial,origin='lower')\n",
+ "ax[1].imshow(gammat_radial,origin='lower')\n",
+ "\n",
+ "ax[0].set_ylabel('r')\n",
+ "ax[0].set_xlabel('$\\phi$')\n",
+ "\n",
+ "ax[1].set_ylabel('r')\n",
+ "ax[1].set_xlabel('$\\phi$')\n",
+ "\n",
+ "ax[0].set_title(r'$\\kappa$')\n",
+ "ax[1].set_title(r'$\\gamma_t$')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1086e4de",
+ "metadata": {},
+ "source": [
+ "We can recover the 1D radial convergence and tangential shear profiles as follows"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 260,
+ "id": "5c6174be",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "kappa_1D_radial = kappa_radial.mean(axis=-1)\n",
+ "gammat_1D_radial = gammat_radial.mean(axis=-1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ebb41586",
+ "metadata": {},
+ "source": [
+ "Plot the 1D radial profiles"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 240,
+ "id": "29a26710",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "r_bins_mid = 0.5 * (r_bins[1:] + r_bins[:-1])\n",
+ "plt.plot(r_bins_mid,kappa_1D_radial,label='$\\\\kappa(r)$')\n",
+ "plt.plot(r_bins_mid,gammat_1D_radial,label='$\\gamma_t(r)$')\n",
+ "plt.legend()\n",
+ "plt.xlabel('r [Mpc]')\n",
+ "plt.ylabel('Lensing amplitude')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "085bd120-6283-4d8f-adeb-bf4e2c781ba2",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "id": "81e0ceb4",
+ "metadata": {},
+ "source": [
+ "## Miscentering\n",
+ "\n",
+ "We now repeat the above procedure, but this time we use an incorrect center (relative to the true hale center) when performing the calculations, the exemplify the impact of miscentering"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 241,
+ "id": "13660944",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n"
+ ]
+ }
+ ],
+ "source": [
+ "# First we define a miscentring parameter\n",
+ "R_mis = 0.05 # Mpc\n",
+ "\n",
+ "# find new center based on R_mis\n",
+ "_phi_Rmis = 2*np.pi*np.random.random() # isotropic angle for mis-centering\n",
+ "\n",
+ "# mis-centered coordinates\n",
+ "center_1 = R_mis*np.cos(_phi_Rmis)\n",
+ "center_2 = R_mis*np.sin(_phi_Rmis)\n",
+ "\n",
+ "# get tangential and radial profiles\n",
+ "et_map, ex_map, radius_map, angle_map = mod.get_Tangential_and_Cross(e1_map, e2_map, [center_1, center_2], dx=dx.value)\n",
+ "kappa_radial, gammat_radial = mod.get_Radial(radius_map,angle_map,\n",
+ " r_bins.value,phi_bins,\n",
+ " kappa_map,et_map)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 305,
+ "id": "1cb19534",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_maps_and_profiles(kappa_map,Sigma_map,e1_map,e2_map,et_map,ex_map,\n",
+ " kappa_radial,gammat_radial,r_bins):\n",
+ " '''quickly plot the convergence, e1,e2, ex,et and radial profiles\n",
+ " '''\n",
+ " fig, ax1 = plt.subplots(3,2,figsize=(10,18))\n",
+ " ax = ax1.flatten()\n",
+ "\n",
+ " plt.sca(ax[0])\n",
+ " plt.imshow(kappa_map,origin = 'lower')\n",
+ " plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ " plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ " plt.colorbar()\n",
+ " plt.title(r'$\\kappa$')\n",
+ "\n",
+ " plt.sca(ax[1])\n",
+ " plt.imshow(Sigma_map,origin = 'lower')\n",
+ " plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ " plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ " plt.colorbar()\n",
+ " plt.title(r'$\\Sigma [M_\\odot\\; Mpc^{-2}]$')\n",
+ "\n",
+ "\n",
+ " plt.sca(ax[2])\n",
+ " plt.imshow(e1_map)\n",
+ " plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ " plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ " plt.colorbar()\n",
+ " plt.title('$e_1$')\n",
+ " \n",
+ " plt.sca(ax[3])\n",
+ " plt.imshow(e2_map)\n",
+ " plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ " plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ " plt.colorbar()\n",
+ " plt.title('$e_2$')\n",
+ " \n",
+ " \n",
+ " plt.sca(ax[4])\n",
+ " plt.imshow(et_map)\n",
+ " plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ " plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ " plt.colorbar()\n",
+ " plt.title('$e_t$')\n",
+ " \n",
+ " plt.sca(ax[5])\n",
+ " plt.imshow(ex_map)\n",
+ " plt.xlabel('x ({0})'.format(r_max.unit))\n",
+ " plt.ylabel('y ({0})'.format(r_max.unit))\n",
+ " plt.colorbar()\n",
+ " plt.title('$e_x$')\n",
+ "\n",
+ " plt.subplots_adjust(bottom=0.3)\n",
+ " plt.show()\n",
+ "\n",
+ "\n",
+ " r_bins_mid = 0.5 * (r_bins[1:] + r_bins[:-1])\n",
+ "\n",
+ " fig,ax = plt.subplots(1,2,figsize=(10,4))\n",
+ " \n",
+ " ax[0].imshow(kappa_radial,origin='lower')\n",
+ " ax[1].imshow(gammat_radial,origin='lower')\n",
+ " \n",
+ " ax[0].set_ylabel('r')\n",
+ " ax[0].set_xlabel('$\\phi$')\n",
+ " \n",
+ " ax[1].set_ylabel('r')\n",
+ " ax[1].set_xlabel('$\\phi$')\n",
+ " \n",
+ " ax[0].set_title(r'$\\kappa$')\n",
+ " ax[1].set_title(r'$\\gamma_t$') \n",
+ " plt.show()\n",
+ "\n",
+ "\n",
+ " plt.figure()\n",
+ " r_bins_mid = 0.5 * (r_bins[1:] + r_bins[:-1])\n",
+ " plt.plot(r_bins_mid,kappa_1D_radial,label='$\\\\kappa(r)$')\n",
+ " plt.plot(r_bins_mid,gammat_1D_radial,label='$\\gamma_t(r)$')\n",
+ " plt.legend()\n",
+ " plt.xlabel('r [Mpc]')\n",
+ " plt.ylabel('Lensing amplitude')\n",
+ "\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e9720a7d",
+ "metadata": {},
+ "source": [
+ "we can plot all of the results at once using the plot_maps_and_profiles function"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 306,
+ "id": "32772d61",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot_maps_and_profiles(kappa_map,Sigma_map,e1_map,e2_map,et_map,ex_map,\n",
+ " kappa_radial,gammat_radial,r_bins) "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "510b4c9c",
+ "metadata": {},
+ "source": [
+ "We can also investigate how the 1D radial lensing profiles change as a function of miscentering"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "53fa70dd",
+ "metadata": {},
+ "source": [
+ "Now define a miscentering distribution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 244,
+ "id": "38b199f1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# define miscentring distribution\n",
+ "R_mis_max = 0.5\n",
+ "N_R_mis = 15\n",
+ "R_mis_arr = np.insert(np.logspace(-2, np.log10(R_mis_max), N_R_mis-1), 0, 0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 245,
+ "id": "c17883c2",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:273: RuntimeWarning: invalid value encountered in divide\n",
+ " tangential shear map\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:324: RuntimeWarning: invalid value encountered in divide\n",
+ "/Users/alessandrafumagalli/DESC/CLMM/clmm/model_profile.py:325: RuntimeWarning: invalid value encountered in divide\n"
+ ]
+ }
+ ],
+ "source": [
+ "et_map = np.zeros((N_R_mis,kappa_map.shape[0],kappa_map.shape[1]))\n",
+ "ex_map = np.zeros(et_map.shape)\n",
+ "radius_map = np.zeros(et_map.shape)\n",
+ "angle_map = np.zeros(et_map.shape)\n",
+ "\n",
+ "kappa_radial = np.zeros((N_R_mis,N_r_bins,N_phi_bins))\n",
+ "gammat_radial = np.zeros(kappa_radial.shape)\n",
+ "\n",
+ "for j in range(len(R_mis_arr)):\n",
+ "\n",
+ " _phi_Rmis = 2*np.pi*np.random.random() # isotropic angle for mis-centering\n",
+ "\n",
+ " # mis-centered coordinates\n",
+ " center_1 = R_mis_arr[j]*np.cos(_phi_Rmis)\n",
+ " center_2 = R_mis_arr[j]*np.sin(_phi_Rmis)\n",
+ "\n",
+ " et_map[j], ex_map[j], radius_map[j], angle_map[j] = mod.get_Tangential_and_Cross(e1_map, e2_map, [center_1, center_2], dx=dx.value)\n",
+ "\n",
+ " kappa_radial[j], gammat_radial[j] = mod.get_Radial(radius_map[j],angle_map[j],\n",
+ " r_bins.value,phi_bins,\n",
+ " kappa_map,et_map[j])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 246,
+ "id": "22ba0469",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "R_mis_ind = 10\n",
+ "plot_maps_and_profiles(kappa_map,e1_map,e2_map,et_map[R_mis_ind],ex_map[R_mis_ind],\n",
+ " kappa_radial[R_mis_ind],gammat_radial[R_mis_ind],r_bins)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2dbde4d3",
+ "metadata": {},
+ "source": [
+ "Finally we can plot the lensing profiles for the whole miscentering distrubtion"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 247,
+ "id": "7e458c10",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "for i in range(N_R_mis):\n",
+ " plt.plot(r_bins_mid,kappa_radial[i].mean(axis=-1),label='misc=%1.2f' %R_mis_arr[i])\n",
+ "plt.xlabel('r [Mpc]')\n",
+ "plt.ylabel('Convergence radial profile')\n",
+ "plt.legend(ncol=3)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 248,
+ "id": "2864956a",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "for i in range(N_R_mis):\n",
+ " plt.plot(r_bins_mid,gammat_radial[i].mean(axis=-1),label='misc=%1.2f' %R_mis_arr[i])\n",
+ "plt.xlabel('r [Mpc]')\n",
+ "plt.ylabel('Gamma_t radial profile')\n",
+ "plt.legend(ncol=3)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ea9ce0dc",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "032de18c-3014-4bdb-a99f-713a22857d33",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3ae1b92d-bc33-4c4d-99b4-8266f38c4709",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}