Repository: LSSTDESC/CLMM (language: python)

Issue:
Allowing different type of `z_src`
In PR #467 we decided it would be good to allow different option for the parameter `z_src`, based on `z_src_model` : float, array or function. 


Discussion / hints:

Fix patch:
diff --git a/clmm/__init__.py b/clmm/__init__.py
index 6e4e83bf1..b2a8dda55 100644
--- a/clmm/__init__.py
+++ b/clmm/__init__.py
@@ -15,4 +15,4 @@
 )
 from . import support
 
-__version__ = '1.4.9'
+__version__ = '1.5.0'
diff --git a/clmm/support/mock_data.py b/clmm/support/mock_data.py
index 09f53e683..e56732eab 100644
--- a/clmm/support/mock_data.py
+++ b/clmm/support/mock_data.py
@@ -7,7 +7,8 @@
 from ..gcdata import GCData
 from ..theory import compute_tangential_shear, compute_convergence
 from ..utils import (convert_units, compute_lensed_ellipticity, validate_argument,
-                     _chang_z_distrib, _srd_z_distrib, _draw_random_points_from_distribution)
+                     _draw_random_points_from_distribution)
+from .. import z_distributions as zdist
 
 def generate_galaxy_catalog(cluster_m, cluster_z, cluster_c, cosmo, zsrc,
                             cluster_ra=0., cluster_dec=0.,
@@ -227,7 +228,7 @@ def _compute_ngals(ngal_density, field_size, cosmo, cluster_z, zsrc, zsrc_min=No
     if isinstance(zsrc, float):
         ngals = int(ngals)
     elif zsrc in ('chang13', 'desc_srd'):
-        z_distrib_func = _chang_z_distrib if zsrc == 'chang13' else _srd_z_distrib
+        z_distrib_func = zdist.chang2013 if zsrc == 'chang13' else zdist.desc_srd
         # Compute the normalisation for the redshift distribution function (z=[0, inf))
         # z_distrib_func(0, is_cdf=True)=0
         norm = z_distrib_func(np.inf, is_cdf=True)
@@ -235,7 +236,7 @@ def _compute_ngals(ngal_density, field_size, cosmo, cluster_z, zsrc, zsrc_min=No
         prob = (z_distrib_func(zsrc_max, is_cdf=True) - z_distrib_func(zsrc_min, is_cdf=True))/norm
         ngals = int(ngals*prob)
     else:
-        raise ValueError(f"zsrc (={zsrc}) must be float, 'chang13' or 'desc_src'")
+        raise ValueError(f"zsrc (={zsrc}) must be float, 'chang13' or 'desc_srd'")
     return ngals
 
 
@@ -270,7 +271,7 @@ def _generate_galaxy_catalog(cluster_m, cluster_z, cluster_c, cosmo, ngals,
                                     cosmo=cosmo, delta_mdef=delta_so,
                                     halo_profile_model=halo_profile_model,
                                     massdef=massdef,
-                                    z_src_model='single_plane')
+                                    z_src_info='discrete')
 
     gamx = np.zeros(ngals)
 
@@ -279,7 +280,7 @@ def _generate_galaxy_catalog(cluster_m, cluster_z, cluster_c, cosmo, ngals,
                                 z_source=galaxy_catalog['ztrue'], cosmo=cosmo,
                                 delta_mdef=delta_so, halo_profile_model=halo_profile_model,
                                 massdef=massdef,
-                                z_src_model='single_plane')
+                                z_src_info='discrete')
 
     c_cl = SkyCoord(cluster_ra*u.deg, cluster_dec*u.deg, frame='icrs')
     c_gal = SkyCoord(galaxy_catalog['ra']*u.deg, galaxy_catalog['dec']*u.deg, frame='icrs')
@@ -362,12 +363,12 @@ def _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals):
     # Draw zsrc from Chang et al. 2013
     elif zsrc == 'chang13':
         zsrc_list = _draw_random_points_from_distribution(
-            zsrc_min, zsrc_max, ngals, _chang_z_distrib)
+            zsrc_min, zsrc_max, ngals, zdist.chang2013)
 
     # Draw zsrc from the distribution used in the DESC SRD (arxiv:1809.01669)
     elif zsrc == 'desc_srd':
         zsrc_list = _draw_random_points_from_distribution(
-            zsrc_min, zsrc_max, ngals, _srd_z_distrib)
+            zsrc_min, zsrc_max, ngals, zdist.desc_srd)
 
     # Draw zsrc from a uniform distribution between zmin and zmax
     elif zsrc == 'uniform':
diff --git a/clmm/theory/func_layer.py b/clmm/theory/func_layer.py
index c7d37518a..43dc2a75a 100644
--- a/clmm/theory/func_layer.py
+++ b/clmm/theory/func_layer.py
@@ -15,7 +15,7 @@
                        convert_profile_mass_concentration)
 
 __all__ = generic.__all__+['compute_3d_density', 'compute_surface_density',
-                           'compute_excess_surface_density','compute_excess_surface_density_2h', 
+                           'compute_excess_surface_density','compute_excess_surface_density_2h',
                            'compute_surface_density_2h',
                            'compute_critical_surface_density',
                            'compute_tangential_shear', 'compute_convergence',
@@ -23,9 +23,9 @@
                            'compute_magnification_bias']
 
 
-def compute_3d_density(
-        r3d, mdelta, cdelta, z_cl, cosmo, delta_mdef=200,
-        halo_profile_model='nfw', massdef='mean', alpha_ein=None, verbose=False, validate_input=True):
+def compute_3d_density(r3d, mdelta, cdelta, z_cl, cosmo, delta_mdef=200,
+                       halo_profile_model='nfw', massdef='mean', alpha_ein=None,
+                       verbose=False, validate_input=True):
     r"""Retrieve the 3d density :math:`\rho(r)`.
 
     Profiles implemented so far are:
@@ -50,28 +50,29 @@ def compute_3d_density(
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical`;
-            * `virial`;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
 
     alpha_ein : float, optional
-        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only available
-        for the NumCosmo backend
+        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
+        available for the NumCosmo backend
     verbose : boolean, optional
-        If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL backends.
+        If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL
+        backends.
     validate_input : bool, optional
         If True (default), the types of the arguments are checked before proceeding.
 
     Returns
     -------
-    rho : array_like, float
+    rho : numpy.ndarray, float
         3-dimensional mass density in units of :math:`M_\odot\ Mpc^{-3}`
 
     Notes
@@ -95,7 +96,8 @@ def compute_3d_density(
 
 
 def compute_surface_density(r_proj, mdelta, cdelta, z_cl, cosmo, delta_mdef=200,
-                            halo_profile_model='nfw', massdef='mean', alpha_ein=None, verbose=False, validate_input=True):
+                            halo_profile_model='nfw', massdef='mean', alpha_ein=None,
+                            verbose=False, validate_input=True):
     r""" Computes the surface mass density
 
     .. math::
@@ -120,29 +122,30 @@ def compute_surface_density(r_proj, mdelta, cdelta, z_cl, cosmo, delta_mdef=200,
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical` ;
-            * `virial` ;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
 
     alpha_ein : float, optional
-        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only available
-        for the NumCosmo backend
+        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
+        available for the NumCosmo backend
     verbose : boolean, optional
-        If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL backends.
+        If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL
+        backends.
     validate_input : bool, optional
         If True (default), the types of the arguments are checked before proceeding.
 
 
     Returns
     -------
-    sigma : array_like, float
+    sigma : numpy.ndarray, float
         2D projected surface density in units of :math:`M_\odot\ Mpc^{-2}`
 
     Notes
@@ -165,7 +168,8 @@ def compute_surface_density(r_proj, mdelta, cdelta, z_cl, cosmo, delta_mdef=200,
     return sigma
 
 def compute_excess_surface_density(r_proj, mdelta, cdelta, z_cl, cosmo, delta_mdef=200,
-                                   halo_profile_model='nfw', massdef='mean', alpha_ein=None, verbose=False, validate_input=True):
+                                   halo_profile_model='nfw', massdef='mean', alpha_ein=None,
+                                   verbose=False, validate_input=True):
     r""" Computes the excess surface density
 
     .. math::
@@ -193,28 +197,29 @@ def compute_excess_surface_density(r_proj, mdelta, cdelta, z_cl, cosmo, delta_md
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical` - not in cluster_toolkit;
-            * `virial` - not in cluster_toolkit;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
 
     alpha_ein : float, optional
-        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only available
-        for the NumCosmo backend
+        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
+        available for the NumCosmo backend
     verbose : boolean, optional
-        If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL backends.
+        If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL
+        backends.
     validate_input : bool, optional
         If True (default), the types of the arguments are checked before proceeding.
 
     Returns
     -------
-    deltasigma : array_like, float
+    deltasigma : numpy.ndarray, float
         Excess surface density in units of :math:`M_\odot\ Mpc^{-2}`.
     """
     gcm.validate_input = validate_input
@@ -231,11 +236,13 @@ def compute_excess_surface_density(r_proj, mdelta, cdelta, z_cl, cosmo, delta_md
     gcm.validate_input = True
     return deltasigma
 
-def compute_excess_surface_density_2h(r_proj, z_cl, cosmo, halobias=1., lsteps=500, validate_input=True):
+def compute_excess_surface_density_2h(r_proj, z_cl, cosmo, halobias=1., lsteps=500,
+                                      validate_input=True):
     r""" Computes the 2-halo term excess surface density from eq.(13) of Oguri & Hamana (2011)
 
     .. math::
-        \Delta\Sigma_{\rm 2h}(R) = \frac{\rho_m(z)b(M)}{(1 + z)^3D_A(z)^2} \int\frac{ldl}{(2\pi)} P_{\rm mm}(k_l, z)J_2(l\theta)
+        \Delta\Sigma_{\rm 2h}(R) = \frac{\rho_m(z)b(M)}{(1 + z)^3D_A(z)^2} \int\frac{ldl}{(2\pi)}
+        P_{\rm mm}(k_l, z)J_2(l\theta)
 
     where
 
@@ -255,19 +262,20 @@ def compute_excess_surface_density_2h(r_proj, z_cl, cosmo, halobias=1., lsteps=5
     halobias : float, optional
         Value of the halo bias
     lsteps : int, optional
-        Steps for the numerical integration 
+        Steps for the numerical integration
     validate_input: bool
         Validade each input argument
 
     Returns
     -------
-    deltasigma_2h : array_like, float
+    deltasigma_2h : numpy.ndarray, float
         2-halo term excess surface density in units of :math:`M_\odot\ Mpc^{-2}`.
     """
     gcm.validate_input = validate_input
     gcm.set_cosmo(cosmo)
 
-    deltasigma_2h = gcm.eval_excess_surface_density_2h(r_proj, z_cl, halobias=halobias, lsteps=lsteps)
+    deltasigma_2h = gcm.eval_excess_surface_density_2h(r_proj, z_cl,
+                                                       halobias=halobias, lsteps=lsteps)
 
     gcm.validate_input = True
     return deltasigma_2h
@@ -276,7 +284,8 @@ def compute_surface_density_2h(r_proj, z_cl, cosmo, halobias=1, lsteps=500, vali
     r""" Computes the 2-halo term surface density from eq.(13) of Oguri & Hamana (2011)
 
     .. math::
-        \Sigma_{\rm 2h}(R) = \frac{\rho_m(z)b(M)}{(1 + z)^3D_A(z)^2} \int\frac{ldl}{(2\pi)} P_{\rm mm}(k_l, z)J_0(l\theta)
+        \Sigma_{\rm 2h}(R) = \frac{\rho_m(z)b(M)}{(1 + z)^3D_A(z)^2} \int\frac{ldl}{(2\pi)}
+        P_{\rm mm}(k_l, z)J_0(l\theta)
 
     where
 
@@ -296,13 +305,13 @@ def compute_surface_density_2h(r_proj, z_cl, cosmo, halobias=1, lsteps=500, vali
     halobias : float, optional
         Value of the halo bias
     lsteps : int, optional
-        Steps for the numerical integration 
+        Steps for the numerical integration
     validate_input: bool
         Validade each input argument
 
     Returns
     -------
-    sigma_2h : array_like, float
+    sigma_2h : numpy.ndarray, float
         2-halo term surface density in units of :math:`M_\odot\ Mpc^{-2}`.
     """
     gcm.validate_input = validate_input
@@ -313,8 +322,10 @@ def compute_surface_density_2h(r_proj, z_cl, cosmo, halobias=1, lsteps=500, vali
     gcm.validate_input = True
     return sigma_2h
 
-def compute_critical_surface_density(cosmo, z_cluster, z_source=None, use_pdz=False, pzbins=None, pzpdf=None, validate_input=True):
-    r"""Computes either 
+def compute_critical_surface_density(cosmo, z_cluster, z_source=None,
+                                     use_pdz=False, pzbins=None, pzpdf=None,
+                                     validate_input=True):
+    r"""Computes either
 
     - the critical surface density if `use_pdz=False`
 
@@ -324,12 +335,14 @@ def compute_critical_surface_density(cosmo, z_cluster, z_source=None, use_pdz=Fa
     or
 
     - the 'effective critical surface density' if `use_pdz=True`
-     
+
     .. math::
-        \langle \Sigma_{\rm crit}^{-1}\rangle^{-1} = \left(\int \frac{1}{\Sigma_{\rm crit}(z)} p(z) dz\right)^{-1} 
-    
-    where :math:`p(z)` is the source photoz probability density function. 
-    This comes from the maximum likelihood estimator for evaluating a :math:`\Delta\Sigma` profile.
+        \langle \Sigma_{\rm crit}^{-1}\rangle^{-1} = \left(\int \frac{1}{\Sigma_{\rm crit}(z)}
+        p(z) dz\right)^{-1}
+
+    where :math:`p(z)` is the source photoz probability density function.
+    This comes from the maximum likelihood estimator for evaluating a :math:`\Delta\Sigma`
+    profile.
 
 
     Parameters
@@ -341,9 +354,10 @@ def compute_critical_surface_density(cosmo, z_cluster, z_source=None, use_pdz=Fa
     z_source : array_like, float
         Background source galaxy redshift(s)
     use_pdz : bool
-        Flag to use the photoz pdf. If `False` (default), `sigma_c` is computed using the source redshift point estimates `z_source`. 
-        If `True`, `sigma_c` is computed as 1/<1/Sigma_crit>, where the average is performed using the individual galaxy redshift pdf. 
-        In that case, the `pzbins` and `pzpdf` should be specified. 
+        Flag to use the photoz pdf. If `False` (default), `sigma_c` is computed using the source
+        redshift point estimates `z_source`. If `True`, `sigma_c` is computed as 1/<1/Sigma_crit>,
+        where the average is performed using the individual galaxy redshift pdf.
+        In that case, the `pzbins` and `pzpdf` should be specified.
     pzbins : array-like
         Bins where the source redshift pdf is defined
     pzpdf : array-like
@@ -354,19 +368,22 @@ def compute_critical_surface_density(cosmo, z_cluster, z_source=None, use_pdz=Fa
 
     Returns
     -------
-    sigma_c : array_like, float
-        Cosmology-dependent (effective) critical surface density in units of :math:`M_\odot\ Mpc^{-2}`
+    sigma_c : numpy.ndarray, float
+        Cosmology-dependent (effective) critical surface density in units of
+        :math:`M_\odot\ Mpc^{-2}`
     """
 
     gcm.validate_input = validate_input
     gcm.set_cosmo(cosmo)
-    sigma_c = gcm.eval_critical_surface_density(z_cluster, z_src=z_source, use_pdz=use_pdz, pzbins=pzbins, pzpdf=pzpdf)
+    sigma_c = gcm.eval_critical_surface_density(z_cluster, z_src=z_source,
+                                                use_pdz=use_pdz, pzbins=pzbins, pzpdf=pzpdf)
 
     gcm.validate_input = True
     return sigma_c
 
 def compute_tangential_shear(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, delta_mdef=200,
-                             halo_profile_model='nfw', massdef='mean', alpha_ein=None, z_src_model='single_plane', 
+                             halo_profile_model='nfw', massdef='mean', alpha_ein=None,
+                             z_src_info='discrete', beta_kwargs=None,
                              verbose=False, validate_input=True):
     r"""Computes the tangential shear
 
@@ -388,8 +405,9 @@ def compute_tangential_shear(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo,
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float
-        Background source galaxy redshift(s)
+    z_source : array_like, float, function
+        Information on the background source galaxy redshift(s). Value required depends on
+        `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
         CLMM Cosmology object
     delta_mdef : int, optional
@@ -397,67 +415,92 @@ def compute_tangential_shear(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo,
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical` - not in cluster_toolkit;
-            * `virial` - not in cluster_toolkit;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
 
     alpha_ein : float, optional
-        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only available
-        for the NumCosmo backend
-    z_src_model : str, optional
-        Source redshift model, with the following supported options:
-            `single_plane` (default) - all sources at one redshift (if
-            `z_source` is a float) or known individual source galaxy redshifts
-            (if `z_source` is an array and `r_proj` is a float);
+        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
+        availablefor the NumCosmo backend
+    z_src_info : str, optional
+        Type of redshift information provided by the `z_source` argument.
+        The following supported options are:
+
+            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
+              It can be individual redshifts for each source galaxy when `z_source` is an array
+              or all sources are at the same redshift when `z_source` is a float.
+
+            * 'distribution' : A redshift distribution function is provided by `z_source`.
+              `z_source` must be a one dimensional function.
+
+            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
+              `z_source` must be a tuple containing
+              ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+              the lensing efficiency and square of the lensing efficiency averaged over
+              the galaxy redshift distribution repectively.
+
+                .. math::
+                    \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                    {D_{L,\infty}}\right\rangle
+
+                .. math::
+                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                    {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+    beta_kwargs: None, dict
+        Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+        Only used if `z_src_info='distribution'`. Possible keys are:
+
+            * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+              when performing the sum. (default=None)
+            * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+              when performing the sum. (default=10.0)
+            * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+              `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
     verbose : bool, optional
-        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL backends. 
+        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL
+        backends.
     validate_input : bool, optional
         If True (default), the types of the arguments are checked before proceeding.
 
     Returns
     -------
-    gammat : array_like, float
+    gammat : numpy.ndarray, float
         Tangential shear
-
-    Notes
-    -----
-    TODO: Implement `known_z_src` (known individual source galaxy redshifts
-    e.g. discrete case) and `z_src_distribution` (known source redshift
-    distribution e.g. continuous case requiring integration) options for
-    `z_src_model`. We will need :math:`\gamma_\infty` and :math:`\kappa_\infty`
-    for alternative z_src_models using :math:`\beta_s`.
     """
-    if z_src_model == 'single_plane':
-        gcm.validate_input = validate_input
-        gcm.set_cosmo(cosmo)
-        gcm.set_halo_density_profile(
-            halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
-        gcm.set_concentration(cdelta)
-        gcm.set_mass(mdelta)
-        if alpha_ein is not None:
-            gcm.set_einasto_alpha(alpha_ein)
-        if np.min(r_proj) < 1.e-11:
-            raise ValueError(
-                f"Rmin = {np.min(r_proj):.2e} Mpc/h! This value is too small "
-                "and may cause computational issues.")
-
-        gammat = gcm.eval_tangential_shear(r_proj, z_cluster, z_source, verbose=verbose)
-    else:
-        raise ValueError("Unsupported z_src_model")
+
+    gcm.validate_input = validate_input
+    gcm.set_cosmo(cosmo)
+    gcm.set_halo_density_profile(
+        halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
+    gcm.set_concentration(cdelta)
+    gcm.set_mass(mdelta)
+    if alpha_ein is not None:
+        gcm.set_einasto_alpha(alpha_ein)
+    if np.min(r_proj) < 1.e-11:
+        raise ValueError(
+            f"Rmin = {np.min(r_proj):.2e} Mpc/h! This value is too small "
+            "and may cause computational issues.")
+
+    tangential_shear = gcm.eval_tangential_shear(r_proj, z_cluster, z_source,
+                                                 z_src_info=z_src_info, beta_kwargs=beta_kwargs,
+                                                 verbose=verbose)
 
     gcm.validate_input = True
-    return gammat
+    return tangential_shear
 
 
 def compute_convergence(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, delta_mdef=200,
-                        halo_profile_model='nfw', massdef='mean', alpha_ein=None, z_src_model='single_plane',
+                        halo_profile_model='nfw', massdef='mean', alpha_ein=None,
+                        z_src_info='discrete', beta_kwargs=None,
                         verbose=False, validate_input=True):
     r"""Computes the mass convergence
 
@@ -479,8 +522,9 @@ def compute_convergence(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, delt
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float
-        Background source galaxy redshift(s)
+    z_source : array_like, float, function
+        Information on the background source galaxy redshift(s). Value required depends on
+        `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
         CLMM Cosmology object
     delta_mdef : int, optional
@@ -488,81 +532,99 @@ def compute_convergence(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, delt
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical` - not in cluster_toolkit;
-            * `virial` - not in cluster_toolkit;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
 
     alpha_ein : float, optional
-        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only available
-        for the NumCosmo backend
-    z_src_model : str, optional
-        Source redshift model, with the following supported options:
-            `single_plane` (default) - all sources at one redshift (if
-            `z_source` is a float) or known individual source galaxy redshifts
-            (if `z_source` is an array and `r_proj` is a float);
+        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
+        available for the NumCosmo backend
+    z_src_info : str, optional
+        Type of redshift information provided by the `z_source` argument.
+        The following supported options are:
+
+            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
+              It can be individual redshifts for each source galaxy when `z_source` is an array
+              or all sources are at the same redshift when `z_source` is a float.
+
+            * 'distribution' : A redshift distribution function is provided by `z_source`.
+              `z_source` must be a one dimensional function.
+
+            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
+              `z_source` must be a tuple containing
+              ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+              the lensing efficiency and square of the lensing efficiency averaged over
+              the galaxy redshift distribution repectively.
+
+                .. math::
+                    \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                    {D_{L,\infty}}\right\rangle
+
+                .. math::
+                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                    {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+    beta_kwargs: None, dict
+        Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+        Only used if `z_src_info='distribution'`. Possible keys are:
+
+            * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+              when performing the sum. (default=None)
+            * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+              when performing the sum. (default=10.0)
+            * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+              `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
     verbose : bool, optional
-        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL backends. 
+        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL
+        backends.
     validate_input : bool, optional
         If True (default), the types of the arguments are checked before proceeding.
 
     Returns
     -------
-    kappa : array_like, float
+    kappa : numpy.ndarray, float
         Mass convergence, kappa.
 
-    Notes
-    -----
-    TODO: Implement `known_z_src` (known individual source galaxy redshifts
-    e.g. discrete case) and `z_src_distribution` (known source redshift
-    distribution e.g. continuous case requiring integration) options for
-    `z_src_model`. We will need :math:`\gamma_\infty` and :math:`\kappa_\infty`
-    for alternative z_src_models using :math:`\beta_s`.
     """
 
-    if z_src_model == 'single_plane':
-
-        gcm.validate_input = validate_input
-        gcm.set_cosmo(cosmo)
-        gcm.set_halo_density_profile(
-            halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
-        gcm.set_concentration(cdelta)
-        gcm.set_mass(mdelta)
-        if alpha_ein is not None:
-            gcm.set_einasto_alpha(alpha_ein)
-
-        kappa = gcm.eval_convergence(r_proj, z_cluster, z_source, verbose=verbose)
-
-    # elif z_src_model == 'known_z_src': # Discrete case
-    #     raise NotImplementedError('Need to implemnt Beta_s functionality, or average'+\
-    #                               'sigma/sigma_c kappa_t = Beta_s*kappa_inf')
-    # elif z_src_model == 'z_src_distribution': # Continuous ( from a distribution) case
-    #     raise NotImplementedError('Need to implement Beta_s calculation from integrating'+\
-    #                               'distribution of redshifts in each radial bin')
-    else:
-        raise ValueError("Unsupported z_src_model")
-
-    if np.any(np.array(z_source) <= z_cluster):
-        warnings.warn(
-            'Some source redshifts are lower than the cluster redshift.'
-            ' kappa = 0 for those galaxies.')
+    gcm.validate_input = validate_input
+    gcm.set_cosmo(cosmo)
+    gcm.set_halo_density_profile(
+        halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
+    gcm.set_concentration(cdelta)
+    gcm.set_mass(mdelta)
+    if alpha_ein is not None:
+        gcm.set_einasto_alpha(alpha_ein)
+
+    convergence = gcm.eval_convergence(r_proj, z_cluster, z_source, z_src_info=z_src_info,
+                                       beta_kwargs=beta_kwargs, verbose=verbose)
+
+    if z_src_info=='discrete':
+        if np.any(np.array(z_source) <= z_cluster):
+            warnings.warn(
+                'Some source redshifts are lower than the cluster redshift.'
+                ' kappa = 0 for those galaxies.')
 
     gcm.validate_input = True
-    return kappa
+    return convergence
 
 
 def compute_reduced_tangential_shear(
         r_proj, mdelta, cdelta, z_cluster, z_source, cosmo,
-        delta_mdef=200, halo_profile_model='nfw', massdef='mean', alpha_ein=None,
-        z_src_model='single_plane', beta_s_mean=None, beta_s_square_mean=None,
-        z_distrib_func=None, validate_input=True, verbose=False):
-    r"""Computes the reduced tangential shear :math:`g_t = \frac{\gamma_t}{1-\kappa}`.
+        delta_mdef=200, halo_profile_model='nfw', massdef='mean', z_src_info='discrete',
+        approx=None, beta_kwargs=None, alpha_ein=None, validate_input=True, verbose=False):
+    r"""Computes the reduced tangential shear
+
+    .. math::
+        g_t = \frac{\gamma_t}{1-\kappa}
 
     Parameters
     ----------
@@ -574,8 +636,9 @@ def compute_reduced_tangential_shear(
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float
-        Background source galaxy redshift(s)
+    z_source : array_like, float, function
+        Information on the background source galaxy redshift(s). Value required depends on
+        `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
         CLMM Cosmology object
     delta_mdef : int, optional
@@ -583,58 +646,101 @@ def compute_reduced_tangential_shear(
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical` - not in cluster_toolkit;
-            * `virial` - not in cluster_toolkit;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
 
-    alpha_ein : float, optional
-        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only available
-        for the NumCosmo backend
+    z_src_info : str, optional
+        Type of redshift information provided by the `z_source` argument.
+        The following supported options are:
 
-    z_src_model : str, optional
-        Source redshift model, with the following supported options:
+            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
+              It can be individual redshifts for each source galaxy when `z_source` is an array
+              or all sources are at the same redshift when `z_source` is a float.
 
-            * `single_plane` (default): all sources at one redshift (if `z_source` is a float) \
-                or known individual source galaxy redshifts (if `z_source` is an array and \
-                `r_proj` is a float);
-            * `applegate14`: use the equation (6) in Weighing the Giants - III \
-                (Applegate et al. 2014; https://arxiv.org/abs/1208.0605) to evaluate tangential reduced shear;
-            * `schrabback18`: use the equation (12) in Cluster Mass Calibration at High Redshift \
-                (Schrabback et al. 2017; https://arxiv.org/abs/1611.03866) to evaluate tangential reduced shear;
+            * 'distribution' : A redshift distribution function is provided by `z_source`.
+              `z_source` must be a one dimensional function.
 
-        z_distrib_func: one-parameter function
-            Redshift distribution function. This function is used to compute the beta values if they are not provided. The default is the Chang et al (2013) distribution function.
-
-        beta_s_mean: array_like, float, optional
-            Lensing efficiency averaged over the galaxy redshift distribution. If not provided, it will be computed using the default redshift distribution or the one given by the user.
+            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
+              `z_source` must be a tuple containing
+              ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+              the lensing efficiency and square of the lensing efficiency averaged over
+              the galaxy redshift distribution repectively.
 
                 .. math::
-                    \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}{D_{L,\infty}}\right\rangle
-
-        beta_s_square_mean: array_like, float, optional
-            Square of the lensing efficiency averaged over the galaxy redshift distribution. If not provided, it will be computed using the default redshift distribution or the one given by the user.
+                    \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                    {D_{L,\infty}}\right\rangle
 
                 .. math::
-                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}{D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
-    
+                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                    {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+    approx : str, optional
+        Type of computation to be made for reduced tangential shears, options are:
+
+            * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
+              If `z_src_info='discrete'`, full computation is made for each
+              `r_proj, z_source` pair individually. If `z_src_info='distribution'`, reduced
+              tangential shear at each value of `r_proj` is calculated as
+
+              .. math::
+                  g_t
+                  =\left<\frac{\beta_s\gamma_{\infty}}{1-\beta_s\kappa_{\infty}}\right>
+                  =\frac{\int_{z_{min}}^{z_{max}}\frac{\beta_s(z)\gamma_{\infty}}
+                  {1-\beta_s(z)\kappa_{\infty}}N(z)\text{d}z}
+                  {\int_{z_{min}}^{z_{max}} N(z)\text{d}z}
+
+            * 'order1' : Same approach as in Weighing the Giants - III (equation 6 in
+              Applegate et al. 2014; https://arxiv.org/abs/1208.0605). `z_src_info` must be
+              either 'beta', or 'distribution' (that will be used to compute
+              :math:`\langle \beta_s \rangle`)
+
+              .. math::
+                  g_t\approx\frac{\left<\beta_s\right>\gamma_{\infty}}
+                  {1-\left<\beta_s\right>\kappa_{\infty}}
+
+            * 'order2' : Same approach as in Cluster Mass Calibration at High
+              Redshift (equation 12 in Schrabback et al. 2017;
+              https://arxiv.org/abs/1611.03866).
+              `z_src_info` must be either 'beta', or 'distribution' (that will be used
+              to compute :math:`\langle \beta_s \rangle` and :math:`\langle \beta_s^2 \rangle`)
+
+              .. math::
+                  g_t\approx\frac{\left<\beta_s\right>\gamma_{\infty}}
+                  {1-\left<\beta_s\right>\kappa_{\infty}}
+                  \left(1+\left(\frac{\left<\beta_s^2\right>}
+                  {\left<\beta_s\right>^2}-1\right)\left<\beta_s\right>\kappa_{\infty}\right)
+
+    beta_kwargs: None, dict
+        Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+        Only used if `z_src_info='distribution'`. Possible keys are:
+
+            * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+              when performing the sum. (default=None)
+            * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+              when performing the sum. (default=10.0)
+            * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+              `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+    alpha_ein : float, optional
+        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
+        available for the NumCosmo backend
     verbose : bool, optional
-        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL backends. 
+        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL
+        backends.
     validate_input : bool, optional
         If True (default), the types of the arguments are checked before proceeding.
 
-                .. math::
-                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}{D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
-
     Returns
     -------
-    gt : array_like, float
+    gt : numpy.ndarray, float
         Reduced tangential shear
 
     """
@@ -648,17 +754,17 @@ def compute_reduced_tangential_shear(
         gcm.set_einasto_alpha(alpha_ein)
 
     red_tangential_shear = gcm.eval_reduced_tangential_shear(
-        r_proj, z_cluster, z_source, z_src_model, beta_s_mean, beta_s_square_mean, z_distrib_func, verbose=verbose)
+        r_proj, z_cluster, z_source, z_src_info=z_src_info, approx=approx,
+        beta_kwargs=beta_kwargs, verbose=verbose)
 
     gcm.validate_input = True
     return red_tangential_shear
 
-# The magnification is computed taking into account just the tangential shear. This is valid for
-# spherically averaged profiles, e.g., NFW and Einasto (by construction the cross shear is zero).
 
 
 def compute_magnification(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, delta_mdef=200,
-                          halo_profile_model='nfw', massdef='mean', alpha_ein=None, z_src_model='single_plane',
+                          halo_profile_model='nfw', massdef='mean', alpha_ein=None,
+                          z_src_info='discrete', approx=None, beta_kwargs=None,
                           verbose=False, validate_input=True):
     r"""Computes the magnification
 
@@ -675,8 +781,9 @@ def compute_magnification(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, de
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float
-        Background source galaxy redshift(s)
+    z_source : array_like, float, function
+        Information on the background source galaxy redshift(s). Value required depends on
+        `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
         CLMM Cosmology object
     delta_mdef : int, optional
@@ -684,69 +791,123 @@ def compute_magnification(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, de
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical` - not in cluster_toolkit;
-            * `virial` - not in cluster_toolkit;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
 
     alpha_ein : float, optional
-        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only available
-        for the NumCosmo backend
-    z_src_model : str, optional
-        Source redshift model, with the following supported options:
-            `single_plane` (default) - all sources at one redshift (if
-            `z_source` is a float) or known individual source galaxy redshifts
-            (if `z_source` is an array and `r_proj` is a float);
+        If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
+        available for the NumCosmo backend
+    z_src_info : str, optional
+        Type of redshift information provided by the `z_source` argument.
+        The following supported options are:
+
+            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
+              It can be individual redshifts for each source galaxy when `z_source` is an array
+              or all sources are at the same redshift when `z_source` is a float.
+
+            * 'distribution' : A redshift distribution function is provided by `z_source`.
+              `z_source` must be a one dimensional function.
+
+            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
+              `z_source` must be a tuple containing
+              ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+              the lensing efficiency and square of the lensing efficiency averaged over
+              the galaxy redshift distribution repectively.
+
+                .. math::
+                    \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                    {D_{L,\infty}}\right\rangle
+
+                .. math::
+                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                    {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+    approx : str, optional
+        Type of computation to be made for magnifications, options are:
+
+            * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
+              If `z_src_info='discrete'`, full computation is made for each
+              `r_proj, z_source` pair individually. If `z_src_info='distribution'`, magnification
+              at each value of `r_proj` is calculated as
+
+              .. math::
+                  \mu
+                  =\left<\frac{1}{\left(1-\beta_s\kappa_{\infty}\right)^2
+                  -\left(\beta_s\gamma_{\infty}\right)^2}\right>
+                  =\frac{\int_{z_{min}}^{z_{max}}\frac{N(z)\text{d}z}
+                  {\left(1-\beta_s(z)\kappa_{\infty}\right)^2
+                  -\left(\beta_s(z)\gamma_{\infty}\right)^2}}
+                  {\int_{z_{min}}^{z_{max}} N(z)\text{d}z}
+
+            * 'order1' : Uses the weak lensing approximation of the magnification with up to
+              first-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+              `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+              compute :math:`\langle \beta_s \rangle`)
+
+              .. math::
+                  \mu \approx 1 + 2 \left<\beta_s\right>\kappa_{\infty}
+
+            * 'order2' : Uses the weak lensing approximation of the magnification with up to
+              second-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+              `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+              compute :math:`\langle \beta_s \rangle` and :math:`\langle \beta_s^2 \rangle`)
+
+              .. math::
+                  \mu \approx 1 + 2 \left<\beta_s\right>\kappa_{\infty}
+                  + 3 \left<\beta_s^2\right>\kappa_{\infty}^2
+                  + \left<\beta_s^2\right>\gamma_{\infty}^2
+
+    beta_kwargs: None, dict
+        Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+        Only used if `z_src_info='distribution'`. Possible keys are:
+
+            * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+              when performing the sum. (default=None)
+            * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+              when performing the sum. (default=10.0)
+            * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+              `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
     verbose : bool, optional
-        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL backends. 
+        If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and CCL
+        backends.
     validate_input : bool, optional
         If True (default), the types of the arguments are checked before proceeding.
 
     Returns
     -------
-    magnification : array_like, float
+    magnification : numpy.ndarray, float
         Magnification :math:`\mu`.
 
-    Notes
-    -----
-    TODO: Implement `known_z_src` (known individual source galaxy redshifts e.g. discrete case) and
-    `z_src_distribution` (known source redshift distribution e.g. continuous case requiring
-    integration) options for `z_src_model`. We will need :math:`\gamma_\infty` and
-    :math:`\kappa_\infty` for alternative z_src_models using :math:`\beta_s`.
     """
 
-    if z_src_model == 'single_plane':
-
-        gcm.validate_input = validate_input
-        gcm.set_cosmo(cosmo)
-        gcm.set_halo_density_profile(
-            halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
-        gcm.set_concentration(cdelta)
-        gcm.set_mass(mdelta)
-        if alpha_ein is not None:
-            gcm.set_einasto_alpha(alpha_ein)
-
-        magnification = gcm.eval_magnification(r_proj, z_cluster, z_source, verbose=verbose)
-
-    # elif z_src_model == 'known_z_src': # Discrete case
-    #     raise NotImplementedError('Need to implemnt Beta_s functionality, or average'+\
-    #                               'sigma/sigma_c kappa_t = Beta_s*kappa_inf')
-    # elif z_src_model == 'z_src_distribution': # Continuous ( from a distribution) case
-    #     raise NotImplementedError('Need to implement Beta_s calculation from integrating'+\
-    #                               'distribution of redshifts in each radial bin')
-    else:
-        raise ValueError("Unsupported z_src_model")
-
-    if np.any(np.array(z_source) <= z_cluster):
-        warnings.warn(
-            'Some source redshifts are lower than the cluster redshift.'
-            ' magnification = 1 for those galaxies.')
+
+    gcm.validate_input = validate_input
+    gcm.set_cosmo(cosmo)
+    gcm.set_halo_density_profile(
+        halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
+    gcm.set_concentration(cdelta)
+    gcm.set_mass(mdelta)
+    if alpha_ein is not None:
+        gcm.set_einasto_alpha(alpha_ein)
+
+    magnification = gcm.eval_magnification(r_proj, z_cluster, z_source, z_src_info=z_src_info,
+                                           approx=approx, beta_kwargs=beta_kwargs,
+                                           verbose=verbose)
+
+    if z_src_info=='discrete':
+        if np.any(np.array(z_source) <= z_cluster):
+            warnings.warn(
+                'Some source redshifts are lower than the cluster redshift.'
+                ' magnification = 1 for those galaxies.')
 
     gcm.validate_input = True
     return magnification
@@ -755,15 +916,16 @@ def compute_magnification(r_proj, mdelta, cdelta, z_cluster, z_source, cosmo, de
 
 def compute_magnification_bias(r_proj, alpha, mdelta, cdelta, z_cluster, z_source, cosmo,
                                delta_mdef=200, halo_profile_model='nfw', massdef='mean',
-                               z_src_model='single_plane', validate_input=True):
+                               z_src_info='discrete', approx=None, beta_kwargs=None,
+                               verbose=False, validate_input=True):
 
-    r""" Computes magnification bias from magnification :math:`\mu` 
+    r""" Computes magnification bias from magnification :math:`\mu`
     and slope parameter :math:`\alpha` as :
 
     .. math::
         \mu^{\alpha - 1}.
 
-    The alpha parameter depends on the source sample and is computed as the slope of the 
+    The alpha parameter depends on the source sample and is computed as the slope of the
     cummulative numer counts at a given magnitude :
 
     .. math::
@@ -788,8 +950,9 @@ def compute_magnification_bias(r_proj, alpha, mdelta, cdelta, z_cluster, z_sourc
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float
-        Background source galaxy redshift(s)
+    z_source : array_like, float, function
+        Information on the background source galaxy redshift(s). Value required depends on
+        `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
         CLMM Cosmology object
     delta_mdef : int, optional
@@ -797,47 +960,119 @@ def compute_magnification_bias(r_proj, alpha, mdelta, cdelta, z_cluster, z_sourc
     halo_profile_model : str, optional
         Profile model parameterization (letter case independent):
 
-            * `nfw` (default);
-            * `einasto` - valid in numcosmo only;
-            * `hernquist` - valid in numcosmo only;
+            * 'nfw' (default)
+            * 'einasto' - not in cluster_toolkit
+            * 'hernquist' - not in cluster_toolkit
 
     massdef : str, optional
         Profile mass definition, with the following supported options (letter case independent):
 
-            * `mean` (default);
-            * `critical` - not in cluster_toolkit;
-            * `virial` - not in cluster_toolkit;
+            * 'mean' (default)
+            * 'critical'
+            * 'virial'
+
+    z_src_info : str, optional
+        Type of redshift information provided by the `z_source` argument.
+        The following supported options are:
 
-    z_src_model : str, optional
-        Source redshift model, with the following supported options:
-            `single_plane` (default) - all sources at one redshift (if
-            `z_source` is a float) or known individual source galaxy redshifts
-            (if `z_source` is an array and `r_proj` is a float);
+            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
+              It can be individual redshifts for each source galaxy when `z_source` is an array
+              or all sources are at the same redshift when `z_source` is a float.
 
+            * 'distribution' : A redshift distribution function is provided by `z_source`.
+              `z_source` must be a one dimensional function.
+
+            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
+              `z_source` must be a tuple containing
+              ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+              the lensing efficiency and square of the lensing efficiency averaged over
+              the galaxy redshift distribution repectively.
+
+                .. math::
+                    \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                    {D_{L,\infty}}\right\rangle
+
+                .. math::
+                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                    {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+    approx : str, optional
+        Type of computation to be made for magnification biases, options are:
+
+            * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
+              If `z_src_info='discrete'`, full computation is made for each
+              `r_proj, z_source` pair individually. If `z_src_info='distribution'`, magnification
+              bias at each value of `r_proj` is calculated as
+
+              .. math::
+                  \mu^{\alpha-1}
+                  &=\left(\left<\frac{1}{\left(1-\beta_s\kappa_{\infty}\right)^2
+                  -\left(\beta_s\gamma_{\infty}\right)^2}\right>\right)^{\alpha-1}
+                  \\\\
+                  &=\frac{\int_{z_{min}}^{z_{max}}\frac{N(z)\text{d}z}
+                  {\left(\left(1-\beta_s(z)\kappa_{\infty}\right)^2
+                  -\left(\beta_s(z)\gamma_{\infty}\right)^2\right)^{\alpha-1}}}
+                  {\int_{z_{min}}^{z_{max}} N(z)\text{d}z}
+
+            * 'order1' : Uses the weak lensing approximation of the magnification bias with up
+              to first-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+              `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+              compute :math:`\langle \beta_s \rangle`)
+
+              .. math::
+                  \mu^{\alpha-1} \approx
+                  1 + \left(\alpha-1\right)\left(2 \left<\beta_s\right>\kappa_{\infty}\right)
+
+            * 'order2' : Uses the weak lensing approximation of the magnification bias with up
+              to second-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+              `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+              compute :math:`\langle \beta_s \rangle` and :math:`\langle \beta_s^2 \rangle`)
+
+              .. math::
+                  \mu^{\alpha-1} \approx
+                  1 &+ \left(\alpha-1\right)\left(2 \left<\beta_s\right>\kappa_{\infty}\right)
+                  \\\\
+                  &+ \left(\alpha-1\right)\left(\left<\beta_s^2\right>\gamma_{\infty}^2\right)
+                  \\\\
+                  &+ \left(2\alpha-1\right)\left(\alpha-1\right)
+                  \left(\left<\beta_s^2\right>\kappa_{\infty}^2\right)
+
+    beta_kwargs: None, dict
+        Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+        Only used if `z_src_info='distribution'`. Possible keys are:
+
+            * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+              when performing the sum. (default=None)
+            * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+              when performing the sum. (default=10.0)
+            * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+              `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+    validate_input : bool, optional
+        If True (default), the types of the arguments are checked before proceeding.
 
     Returns
     -------
-    magnification_bias : array_like
+    magnification_bias : numpy.ndarray
         magnification bias
     """
-    if np.any(np.array(z_source) <= z_cluster):
-        warnings.warn(
-            'Some source redshifts are lower than the cluster redshift.'
-            ' magnification = 1 for those galaxies.')
-    if z_src_model == 'single_plane':
-
-        gcm.validate_input = validate_input
-        gcm.set_cosmo(cosmo)
-        gcm.set_halo_density_profile(
-            halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
-        gcm.set_concentration(cdelta)
-        gcm.set_mass(mdelta)
 
-        magnification_bias = gcm.eval_magnification_bias(r_proj, z_cluster, z_source, alpha)
+    gcm.validate_input = validate_input
+    gcm.set_cosmo(cosmo)
+    gcm.set_halo_density_profile(
+        halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef)
+    gcm.set_concentration(cdelta)
+    gcm.set_mass(mdelta)
 
-    else:
-        raise ValueError("Unsupported z_src_model")
+    magnification_bias = gcm.eval_magnification_bias(r_proj, z_cluster, z_source, alpha,
+                                                     z_src_info=z_src_info, approx=approx,
+                                                     beta_kwargs=beta_kwargs, verbose=verbose)
 
+    if z_src_info=='discrete':
+        if np.any(np.array(z_source) <= z_cluster):
+            warnings.warn(
+                'Some source redshifts are lower than the cluster redshift.'
+                ' magnification bias = 1 for those galaxies.')
 
     gcm.validate_input = True
     return magnification_bias
diff --git a/clmm/theory/numcosmo.py b/clmm/theory/numcosmo.py
index 27cf5d600..3a7dccfe3 100644
--- a/clmm/theory/numcosmo.py
+++ b/clmm/theory/numcosmo.py
@@ -165,7 +165,7 @@ def _eval_convergence(self, r_proj, z_cl, z_src):
             self.hdpm, self.cosmo.be_cosmo, r_proj, z_src, z_cl, z_cl)
         return np.vectorize(func)(r_proj, z_src, z_cl)
 
-    def _eval_reduced_tangential_shear_sp(self, r_proj, z_cl, z_src):
+    def _eval_reduced_tangential_shear(self, r_proj, z_cl, z_src):
         """"eval reduced tangential shear considering a single redshift plane for background sources"""
 
         self.cosmo.smd.prepare_if_needed(self.cosmo.be_cosmo)
diff --git a/clmm/theory/parent_class.py b/clmm/theory/parent_class.py
index 728d395e2..609205b2e 100644
--- a/clmm/theory/parent_class.py
+++ b/clmm/theory/parent_class.py
@@ -4,7 +4,7 @@
 import numpy as np
 
 # functions for the 2h term
-from scipy.integrate import simps
+from scipy.integrate import simps, quad
 from scipy.special import jv
 from scipy.interpolate import interp1d
 
@@ -12,7 +12,7 @@
                       compute_magnification_bias_from_magnification,
                       compute_rdelta, compute_profile_mass_in_radius,
                       convert_profile_mass_concentration)
-from ..utils import validate_argument, _integ_pzfuncs, compute_beta_s_mean, compute_beta_s_square_mean
+from ..utils import validate_argument, _integ_pzfuncs, compute_beta_s_mean, compute_beta_s_square_mean, compute_beta_s_func
 
 
 class CLMModeling:
@@ -168,19 +168,16 @@ def _eval_excess_surface_density(self, r_proj, z_cl):
     # 3. Functions that can be used by all subclasses
 
 
-    def _eval_critical_surface_density(self, z_len, z_src):
-        return self.cosmo.eval_sigma_crit(z_len, z_src)
-
     def _set_cosmo(self, cosmo):
         r""" Sets the cosmology to the internal cosmology object"""
         self.cosmo = cosmo if cosmo is not None else self.cosmo_class()
 
-    def _eval_excess_surface_density_2h(self, r_proj, z_cl, halobias=1.,lsteps=500):
+    def _eval_excess_surface_density_2h(self, r_proj, z_cl, halobias=1., lsteps=500):
         """"eval excess surface density from the 2-halo term"""
         da = self.cosmo.eval_da(z_cl)
         rho_m = self.cosmo._get_rho_m(z_cl)
 
-        kk = np.logspace(-5.,5.,1000)
+        kk = np.logspace(-5., 5., 1000)
         pk = self.cosmo._eval_linear_matter_powerspectrum(kk, z_cl)
         interp_pk = interp1d(kk, pk, kind='cubic')
         theta = r_proj / da
@@ -199,7 +196,7 @@ def _eval_surface_density_2h(self, r_proj, z_cl, halobias=1., lsteps=500):
         da = self.cosmo.eval_da(z_cl)
         rho_m = self.cosmo._get_rho_m(z_cl)
 
-        kk = np.logspace(-5.,5.,1000)
+        kk = np.logspace(-5., 5., 1000)
         pk = self.cosmo._eval_linear_matter_powerspectrum(kk, z_cl)
         interp_pk = interp1d(kk, pk, kind='cubic')
         theta = r_proj / da
@@ -213,6 +210,32 @@ def __integrand__( l , theta ):
         val = np.array( [ simps( __integrand__( ll , t ) , ll ) for t in theta ] )
         return halobias * val * rho_m / ( 2 * np.pi  * ( 1 + z_cl )**3 * da**2 )
 
+    def _eval_rdelta(self, z_cl):
+        return compute_rdelta(self.mdelta, z_cl, self.cosmo, self.massdef, self.delta_mdef)
+
+    def _eval_mass_in_radius(self, r3d, z_cl):
+        alpha = self._get_einasto_alpha(z_cl) if self.halo_profile_model=='einasto' else None
+        return compute_profile_mass_in_radius(
+            r3d, z_cl, self.cosmo, self.mdelta, self.cdelta,
+            self.massdef, self.delta_mdef, self.halo_profile_model, alpha)
+
+    def _convert_mass_concentration(self, z_cl, massdef=None, delta_mdef=None,
+                                    halo_profile_model=None, alpha=None):
+        alpha1 = self._get_einasto_alpha(z_cl) if self.halo_profile_model=='einasto' else None
+        return convert_profile_mass_concentration(
+            self.mdelta, self.cdelta, z_cl, self.cosmo,
+            massdef=self.massdef, delta_mdef=self.delta_mdef,
+            halo_profile_model=self.halo_profile_model, alpha=alpha1,
+            massdef2=massdef, delta_mdef2=delta_mdef,
+            halo_profile_model2=halo_profile_model, alpha2=alpha)
+
+
+    # 3.1. All these functions are for the single plane case
+
+
+    def _eval_critical_surface_density(self, z_len, z_src):
+        return self.cosmo.eval_sigma_crit(z_len, z_src)
+
     def _eval_tangential_shear(self, r_proj, z_cl, z_src):
         delta_sigma = self.eval_excess_surface_density(r_proj, z_cl)
         sigma_c = self.eval_critical_surface_density(z_cl, z_src)
@@ -223,7 +246,7 @@ def _eval_convergence(self, r_proj, z_cl, z_src, verbose=False):
         sigma_c = self.eval_critical_surface_density(z_cl, z_src)
         return sigma/sigma_c
 
-    def _eval_reduced_tangential_shear_sp(self, r_proj, z_cl, z_src):
+    def _eval_reduced_tangential_shear(self, r_proj, z_cl, z_src):
         kappa = self.eval_convergence(r_proj, z_cl, z_src)
         gamma_t = self.eval_tangential_shear(r_proj, z_cl, z_src)
         return compute_reduced_shear_from_convergence(gamma_t, kappa)
@@ -237,25 +260,6 @@ def _eval_magnification_bias(self, r_proj, z_cl, z_src, alpha):
         magnification = self.eval_magnification(r_proj, z_cl, z_src)
         return compute_magnification_bias_from_magnification(magnification, alpha)
 
-    def _eval_rdelta(self, z_cl):
-        return compute_rdelta(self.mdelta, z_cl, self.cosmo, self.massdef, self.delta_mdef)
-
-    def _eval_mass_in_radius(self, r3d, z_cl):
-        alpha = self._get_einasto_alpha(z_cl) if self.halo_profile_model=='einasto' else None
-        return compute_profile_mass_in_radius(
-            r3d, z_cl, self.cosmo, self.mdelta, self.cdelta,
-            self.massdef, self.delta_mdef, self.halo_profile_model, alpha)
-
-    def _convert_mass_concentration(self, z_cl, massdef=None, delta_mdef=None,
-                                    halo_profile_model=None, alpha=None):
-        alpha1 = self._get_einasto_alpha(z_cl) if self.halo_profile_model=='einasto' else None
-        return convert_profile_mass_concentration(
-            self.mdelta, self.cdelta, z_cl, self.cosmo,
-            massdef=self.massdef, delta_mdef=self.delta_mdef,
-            halo_profile_model=self.halo_profile_model, alpha=alpha1,
-            massdef2=massdef, delta_mdef2=delta_mdef,
-            halo_profile_model2=halo_profile_model, alpha2=alpha)
-
 
     # 4. Wrapper functions for input validation
 
@@ -417,7 +421,7 @@ def eval_critical_surface_density(self, z_len, z_src=None, use_pdz=False, pzbins
         z_src : array_like, float
             Background source galaxy redshift(s)
         use_pdz : bool
-            Flag to use the photoz pdf. If `False` (default), `sigma_c` is computed using the source redshift point estimates `z_source`.
+            Flag to use the photoz pdf. If `False` (default), `sigma_c` is computed using the source redshift point estimates `z_src`.
             If `True`, `sigma_c` is computed as 1/<1/Sigma_crit>, where the average is performed using
             the individual galaxy redshift pdf. In that case, the `pzbins` and `pzpdf` should be specified.
 
@@ -446,7 +450,7 @@ def eval_critical_surface_density(self, z_len, z_src=None, use_pdz=False, pzbins
                 raise ValueError('Redshift bins and source redshift pdf must be provided when use_pdz is True')
             else:
                 def inv_sigmac(redshift):
-                    return 1./self._eval_critical_surface_density(z_len=z_len, z_src=redshift)        
+                    return 1./self._eval_critical_surface_density(z_len=z_len, z_src=redshift)
                 return 1./_integ_pzfuncs(pzpdf, pzbins, kernel=inv_sigmac, is_unique_pzbins=np.all(pzbins==pzbins[0]))
 
     def eval_surface_density(self, r_proj, z_cl, verbose=False):
@@ -587,7 +591,132 @@ def eval_surface_density_2h(self, r_proj, z_cl, halobias=1., lsteps=500):
         else:
             return self._eval_surface_density_2h(r_proj, z_cl, halobias=halobias, lsteps=lsteps)
 
-    def eval_tangential_shear(self, r_proj, z_cl, z_src, verbose=False):
+    def _get_beta_s_mean(self, z_cl, z_src, z_src_info='discrete', beta_kwargs=None):
+        r"""Get mean value of the geometric lensing efficicency ratio from typical class function.
+
+        Parameters
+        ----------
+        z_cl : float
+            Galaxy cluster redshift
+        z_src : array_like, float, function
+            Information on the background source galaxy redshift(s). Value required depends on
+            `z_src_info` (see below).
+        z_src_info : str, optional
+            Type of redshift information provided by the `z_src` argument.
+            The following supported options are:
+
+                * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+                  It can be individual redshifts for each source galaxy when `z_src` is an
+                  arrayor all sources are at the same redshift when `z_src` is a float.
+
+                * 'distribution' : A redshift distribution function is provided by `z_src`.
+                  `z_src` must be a one dimentional function.
+
+                * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+                  `z_src` must be a tuple containing
+                  ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+                  the lensing efficiency and square of the lensing efficiency averaged over
+                  the galaxy redshift distribution repectively.
+
+                    .. math::
+                        \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                        {D_{L,\infty}}\right\rangle
+
+                    .. math::
+                        \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                        {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+        beta_kwargs: None, dict
+            Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+            Only used if `z_src_info='distribution'`. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        Returns
+        -------
+        array_like, float
+            The averaged lensing efficiency.
+        """
+        if z_src_info=='beta':
+            # z_src (tuple) is (beta_s_mean, beta_s_square_mean)
+            beta_s_mean = z_src[0]
+        elif z_src_info=='distribution':
+            # z_src (function) if PDZ
+            beta_kwargs = {} if beta_kwargs is None else beta_kwargs
+            beta_s_mean = compute_beta_s_mean(
+                z_cl, self.z_inf, self.cosmo, z_distrib_func=z_src, **beta_kwargs)
+        return beta_s_mean
+
+    def _get_beta_s_square_mean(self, z_cl, z_src, z_src_info='discrete',
+                                beta_kwargs=None):
+        r"""Get mean value of the square geometric lensing efficicency ratio from typical class
+        function.
+
+        Parameters
+        ----------
+        z_cl : float
+            Galaxy cluster redshift
+        z_src : array_like, float, function
+            Information on the background source galaxy redshift(s). Value required depends on
+            `z_src_info` (see below).
+        z_src_info : str, optional
+            Type of redshift information provided by the `z_src` argument.
+            The following supported options are:
+
+                * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+                  It can be individual redshifts for each source galaxy when `z_src` is an
+                  arrayor all sources are at the same redshift when `z_src` is a float.
+
+                * 'distribution' : A redshift distribution function is provided by `z_src`.
+                  `z_src` must be a one dimentional function.
+
+                * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+                  `z_src` must be a tuple containing
+                  ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+                  the lensing efficiency and square of the lensing efficiency averaged over
+                  the galaxy redshift distribution repectively.
+
+                    .. math::
+                        \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                        {D_{L,\infty}}\right\rangle
+
+                    .. math::
+                        \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                        {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+        beta_kwargs: None, dict
+            Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+            Only used if `z_src_info='distribution'`. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        Returns
+        -------
+        array_like, float
+            The square averaged lensing efficiency.
+        """
+        if z_src_info=='beta':
+            # z_src (tuple) is (beta_s_mean, beta_s_square_mean)
+            beta_s_square_mean = z_src[1]
+        elif z_src_info=='distribution':
+            # z_src (function) if PDZ
+            beta_kwargs = {} if beta_kwargs is None else beta_kwargs
+            beta_s_square_mean = compute_beta_s_square_mean(
+                z_cl, self.z_inf, self.cosmo, z_distrib_func=z_src, **beta_kwargs)
+        return beta_s_square_mean
+
+    def eval_tangential_shear(self, r_proj, z_cl, z_src, z_src_info='discrete',
+                              beta_kwargs=None, verbose=False):
         r"""Computes the tangential shear
 
         Parameters
@@ -596,8 +725,48 @@ def eval_tangential_shear(self, r_proj, z_cl, z_src, verbose=False):
             The projected radial positions in :math:`M\!pc`.
         z_cl : float
             Galaxy cluster redshift
-        z_src : array_like, float
-            Background source galaxy redshift(s)
+        z_src : array_like, float, function
+            Information on the background source galaxy redshift(s). Value required depends on
+            `z_src_info` (see below).
+        z_src_info : str, optional
+            Type of redshift information provided by the `z_src` argument.
+            The following supported options are:
+
+                * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+                  It can be individual redshifts for each source galaxy when `z_src` is an
+                  arrayor all sources are at the same redshift when `z_src` is a float.
+
+                * 'distribution' : A redshift distribution function is provided by `z_src`.
+                  `z_src` must be a one dimentional function.
+
+                * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+                  `z_src` must be a tuple containing
+                  ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+                  the lensing efficiency and square of the lensing efficiency averaged over
+                  the galaxy redshift distribution repectively.
+
+                    .. math::
+                        \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                        {D_{L,\infty}}\right\rangle
+
+                    .. math::
+                        \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                        {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+        beta_kwargs: None, dict
+            Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+            Only used if `z_src_info='distribution'`. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        verbose : bool, optional
+            If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
+            CCL backends.
 
         Returns
         -------
@@ -607,14 +776,29 @@ def eval_tangential_shear(self, r_proj, z_cl, z_src, verbose=False):
         if self.validate_input:
             validate_argument(locals(), 'r_proj', 'float_array', argmin=0)
             validate_argument(locals(), 'z_cl', float, argmin=0)
-            validate_argument(locals(), 'z_src', 'float_array', argmin=0)
+            validate_argument(locals(), 'z_src_info', str)
+            self._validate_z_src(locals())
 
         if self.halo_profile_model=='einasto' and verbose:
             print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
 
-        return self._eval_tangential_shear(r_proj=r_proj, z_cl=z_cl, z_src=z_src)
+        if z_src_info=='discrete':
+            gammat = self._eval_tangential_shear(r_proj=r_proj, z_cl=z_cl, z_src=z_src)
+        elif z_src_info in ('distribution', 'beta'):
+            beta_s_mean = self._get_beta_s_mean(
+                z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
+
+            gammat_inf = self._eval_tangential_shear(r_proj=r_proj, z_cl=z_cl, z_src=self.z_inf)
+
+            gammat = beta_s_mean * gammat_inf
+        else:
+            raise ValueError(f"Unsupported z_src_info (='{z_src_info}')")
+
+        return gammat
+
+    def eval_convergence(self, r_proj, z_cl, z_src, z_src_info='discrete',
+                         beta_kwargs=None, verbose=False):
 
-    def eval_convergence(self, r_proj, z_cl, z_src, verbose=False):
         r"""Computes the mass convergence
 
         .. math::
@@ -623,7 +807,7 @@ def eval_convergence(self, r_proj, z_cl, z_src, verbose=False):
         or
 
         .. math::
-            \kappa = \kappa_\infty \times \beta_s
+            \kappa = \kappa_\infty \times <\beta_s>
 
         Parameters
         ----------
@@ -631,8 +815,48 @@ def eval_convergence(self, r_proj, z_cl, z_src, verbose=False):
             The projected radial positions in :math:`M\!pc`.
         z_cl : float
             Galaxy cluster redshift
-        z_src : array_like, float
-            Background source galaxy redshift(s)
+        z_src : array_like, float, function
+            Information on the background source galaxy redshift(s). Value required depends on
+            `z_src_info` (see below).
+        z_src_info : str, optional
+            Type of redshift information provided by the `z_src` argument.
+            The following supported options are:
+
+                * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+                  It can be individual redshifts for each source galaxy when `z_src` is an
+                  array or all sources are at the same redshift when `z_src` is a float.
+
+                * 'distribution' : A redshift distribution function is provided by `z_src`.
+                  `z_src` must be a one dimentional function.
+
+                * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+                  `z_src` must be a tuple containing
+                  ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+                  the lensing efficiency and square of the lensing efficiency averaged over
+                  the galaxy redshift distribution repectively.
+
+                    .. math::
+                        \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                        {D_{L,\infty}}\right\rangle
+
+                    .. math::
+                        \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                        {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+        beta_kwargs: None, dict
+            Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+            Only used if `z_src_info='distribution'`. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        verbose : bool, optional
+            If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
+            CCL backends.
 
         Returns
         -------
@@ -642,46 +866,161 @@ def eval_convergence(self, r_proj, z_cl, z_src, verbose=False):
         if self.validate_input:
             validate_argument(locals(), 'r_proj', 'float_array', argmin=0)
             validate_argument(locals(), 'z_cl', float, argmin=0)
-            validate_argument(locals(), 'z_src', 'float_array', argmin=0)
+            validate_argument(locals(), 'z_src_info', str)
+            self._validate_z_src(locals())
 
         if self.halo_profile_model=='einasto' and verbose:
             print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
 
-        return self._eval_convergence(r_proj=r_proj, z_cl=z_cl, z_src=z_src)
+        if z_src_info=='discrete':
+            kappa = self._eval_convergence(r_proj=r_proj, z_cl=z_cl, z_src=z_src)
+        elif z_src_info in ('distribution', 'beta'):
+            beta_s_mean = self._get_beta_s_mean(
+                z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
 
-    def eval_reduced_tangential_shear(self, r_proj, z_cl, z_src, z_src_model='single_plane',
-                                      beta_s_mean=None, beta_s_square_mean=None, z_distrib_func=None, verbose=False):
-        r"""Computes the reduced tangential shear :math:`g_t = \frac{\gamma_t}{1-\kappa}`.
+            kappa_inf = self._eval_convergence(r_proj=r_proj, z_cl=z_cl, z_src=self.z_inf)
+
+            kappa = beta_s_mean * kappa_inf
+        else:
+            raise ValueError(f"Unsupported z_src_info (='{z_src_info}')")
+
+        return kappa
+
+    def _pdz_weighted_avg(self, core, pdz_func, r_proj, z_cl, integ_kwargs=None):
+        r"""Computes function averaged over PDZ
 
         Parameters
         ----------
+        core : function
+            Function to be averaged, must take tangential shear and convergence as input.
+        pdz_func : function
+            Redshift distribution function. Must be a one dimentional function.
         r_proj : array_like
             The projected radial positions in :math:`M\!pc`.
         z_cl : float
             Galaxy cluster redshift
-        z_src : array_like, float
-            Background source galaxy redshift(s)
-        z_src_model : str, optional
-            Source redshift model, with the following supported options:
+        integ_kwargs: None, dict
+            Extra arguments for the redshift integration. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        Returns
+        -------
+        array_like
+            Function averaged by pdz, with r_proj dimention.
+        """
+        tfunc = lambda z, r: compute_beta_s_func(
+            z, z_cl, self.z_inf, self.cosmo, self._eval_tangential_shear, r, z_cl, self.z_inf)
+        kfunc = lambda z, r: compute_beta_s_func(
+            z, z_cl, self.z_inf, self.cosmo, self._eval_convergence, r, z_cl, self.z_inf)
+        __integrand__ = lambda z, r: pdz_func(z)*core(tfunc(z, r), kfunc(z, r))
 
-                * `single_plane` (default): all sources at one redshift (if `z_source` is a float)    or known individual source galaxy redshifts (if `z_source` is an array and `r_proj` is a float).
-                * `applegate14`: use the equation (6) in Weighing the Giants - III (Applegate et al. 2014; https://arxiv.org/abs/1208.0605) to evaluate tangential reduced shear.
-                * `schrabback18`: use the equation (12) in Cluster Mass Calibration at High Redshift (Schrabback et al. 2017; https://arxiv.org/abs/1611.03866) to evaluate tangential reduced shear.
+        _integ_kwargs = {'zmax': 10.0, 'delta_z_cut': 0.1}
+        _integ_kwargs.update({} if integ_kwargs is None else integ_kwargs)
 
-        z_distrib_func: one-parameter function
-            Redshift distribution function. This function is used to compute the beta values if they are not provided. The default is the Chang et al (2013) distribution function.
+        zmax = _integ_kwargs['zmax']
+        delta_z_cut = _integ_kwargs['delta_z_cut']
+        zmin = _integ_kwargs.get('zmin', z_cl+delta_z_cut)
 
-        beta_s_mean: array_like, float, optional
-            Lensing efficiency averaged over the galaxy redshift distribution. If not provided, it will be computed using the default redshift distribution or the one given by the user.
+        out = np.array([quad(__integrand__, zmin, zmax, (r))[0] for r in r_proj])
+        return out/quad(pdz_func, zmin, zmax)[0]
 
-                .. math::
-                    \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}{D_{L,\infty}}\right\rangle
+    def eval_reduced_tangential_shear(self, r_proj, z_cl, z_src, z_src_info='discrete',
+                                      approx=None, beta_kwargs=None, verbose=False):
+        r"""Computes the reduced tangential shear
 
-        beta_s_square_mean: array_like, float, optional
-            Square of the lensing efficiency averaged over the galaxy redshift distribution. If not provided, it will be computed using the default redshift distribution or the one given by the user.
+        .. math::
+            g_t = \frac{\gamma_t}{1-\kappa}
 
-                .. math::
-                    \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}{D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+        Parameters
+        ----------
+        r_proj : array_like
+            The projected radial positions in :math:`M\!pc`.
+        z_cl : float
+            Galaxy cluster redshift
+        z_src : array_like, float, function
+            Information on the background source galaxy redshift(s). Value required depends on
+            `z_src_info` (see below).
+        z_src_info : str, optional
+            Type of redshift information provided by the `z_src` argument.
+            The following supported options are:
+
+                * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+                  It can be individual redshifts for each source galaxy when `z_src` is an
+                  array or all sources are at the same redshift when `z_src` is a float.
+
+                * 'distribution' : A redshift distribution function is provided by `z_src`.
+                  `z_src` must be a one dimentional function.
+
+                * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+                  `z_src` must be a tuple containing
+                  ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+                  the lensing efficiency and square of the lensing efficiency averaged over
+                  the galaxy redshift distribution repectively.
+
+                    .. math::
+                        \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                        {D_{L,\infty}}\right\rangle
+
+                    .. math::
+                        \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                        {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+        approx : str, optional
+            Type of computation to be made for reduced tangential shears, options are:
+
+                * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
+                  If `z_src_info='discrete'`, full computation is made for each
+                  `r_proj, z_src` pair individually. If `z_src_info='distribution'`, reduced
+                  tangential shear at each value of `r_proj` is calculated as
+
+                  .. math::
+                      g_t
+                      =\left<\frac{\beta_s\gamma_{\infty}}{1-\beta_s\kappa_{\infty}}\right>
+                      =\frac{\int_{z_{min}}^{z_{max}}\frac{\beta_s(z)\gamma_{\infty}}
+                      {1-\beta_s(z)\kappa_{\infty}}N(z)\text{d}z}
+                      {\int_{z_{min}}^{z_{max}} N(z)\text{d}z}
+
+                * 'order1' : Same approach as in Weighing the Giants - III (equation 6 in
+                  Applegate et al. 2014; https://arxiv.org/abs/1208.0605). `z_src_info` must be
+                  either 'beta', or 'distribution' (that will be used to compute
+                  :math:`\langle \beta_s \rangle`)
+
+                  .. math::
+                      g_t\approx\frac{\left<\beta_s\right>\gamma_{\infty}}
+                      {1-\left<\beta_s\right>\kappa_{\infty}}
+
+                * 'order2' : Same approach as in Cluster Mass Calibration at High
+                  Redshift (equation 12 in Schrabback et al. 2017;
+                  https://arxiv.org/abs/1611.03866).
+                  `z_src_info` must be either 'beta', or 'distribution' (that will be used
+                  to compute :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`)
+
+                  .. math::
+                      g_t\approx\frac{\left<\beta_s\right>\gamma_{\infty}}
+                      {1-\left<\beta_s\right>\kappa_{\infty}}
+                      \left(1+\left(\frac{\left<\beta_s^2\right>}
+                      {\left<\beta_s\right>^2}-1\right)\left<\beta_s\right>\kappa_{\infty}\right)
+
+        beta_kwargs: None, dict
+            Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+            Only used if `z_src_info='distribution'`. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        verbose : bool, optional
+            If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
+            CCL backends.
 
         Returns
         -------
@@ -695,35 +1034,47 @@ def eval_reduced_tangential_shear(self, r_proj, z_cl, z_src, z_src_model='single
         if self.validate_input:
             validate_argument(locals(), 'r_proj', 'float_array', argmin=0)
             validate_argument(locals(), 'z_cl', float, argmin=0)
-            validate_argument(locals(), 'z_src', 'float_array', argmin=0)
+            validate_argument(locals(), 'z_src_info', str)
+            validate_argument(locals(), 'approx', str, none_ok=True)
+            self._validate_z_src(locals())
 
         if self.halo_profile_model=='einasto' and verbose:
             print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
 
-        if z_src_model == 'single_plane':
-            gt = self._eval_reduced_tangential_shear_sp(r_proj, z_cl, z_src)
-
-        elif z_src_model == 'applegate14':
-            if beta_s_mean is None or beta_s_square_mean is None:
-                beta_s_mean = compute_beta_s_mean (z_cl, self.z_inf, self.cosmo, z_distrib_func=z_distrib_func)
-                beta_s_square_mean = compute_beta_s_square_mean (z_cl, self.z_inf, self.cosmo, z_distrib_func=z_distrib_func)
-            gammat_inf = self._eval_tangential_shear(r_proj, z_cl, self.z_inf)
-            kappa_inf = self._eval_convergence(r_proj, z_cl, self.z_inf)
-            gt = beta_s_mean * gammat_inf / (1. - beta_s_square_mean / beta_s_mean * kappa_inf)
-
-        elif z_src_model == 'schrabback18':
-            if beta_s_mean is None or beta_s_square_mean is None:
-                beta_s_mean = compute_beta_s_mean (z_cl, self.z_inf, self.cosmo, z_distrib_func=z_distrib_func)
-                beta_s_square_mean = compute_beta_s_square_mean (z_cl, self.z_inf, self.cosmo, z_distrib_func=z_distrib_func)
-            gammat_inf = self._eval_tangential_shear(r_proj, z_cl, self.z_inf)
-            kappa_inf = self._eval_convergence(r_proj, z_cl, self.z_inf)
-            gt = (1. + (beta_s_square_mean / (beta_s_mean * beta_s_mean) - 1.) * beta_s_mean * kappa_inf) * (beta_s_mean * gammat_inf / (1. - beta_s_mean * kappa_inf))
+        if approx is None:
+            if z_src_info=='distribution':
+                core = lambda gammat, kappa: gammat/(1-kappa)
+                gt = self._pdz_weighted_avg(core, z_src, r_proj, z_cl,
+                                            integ_kwargs=beta_kwargs)
+            elif z_src_info=='discrete':
+                gt = self._eval_reduced_tangential_shear(r_proj, z_cl, z_src)
+            else:
+                raise ValueError(
+                    "approx=None requires z_src_info='discrete' or 'distribution',"
+                    f"z_src_info='{z_src_info}' was provided.")
+
+        elif approx in ('order1', 'order2'):
+            beta_s_mean = self._get_beta_s_mean(
+                z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
+
+            gammat_inf = self._eval_tangential_shear(r_proj, z_cl, z_src=self.z_inf)
+            kappa_inf = self._eval_convergence(r_proj, z_cl, z_src=self.z_inf)
+
+            gt = beta_s_mean * gammat_inf / (1. - beta_s_mean * kappa_inf)
 
+            if approx == 'order2':
+                beta_s_square_mean = self._get_beta_s_square_mean(
+                    z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
+
+                gt *= (1. + (beta_s_square_mean / (beta_s_mean * beta_s_mean) - 1.) \
+                           * beta_s_mean * kappa_inf )
         else:
-            raise ValueError("Unsupported z_src_model")
+            raise ValueError(f"Unsupported approx (='{approx}')")
+
         return gt
 
-    def eval_magnification(self, r_proj, z_cl, z_src, verbose=False):
+    def eval_magnification(self, r_proj, z_cl, z_src, z_src_info='discrete',
+                           approx=None, beta_kwargs=None, verbose=False):
         r"""Computes the magnification
 
         .. math::
@@ -735,31 +1086,133 @@ def eval_magnification(self, r_proj, z_cl, z_src, verbose=False):
             The projected radial positions in :math:`M\!pc`.
         z_cl : float
             Galaxy cluster redshift
-        z_src : array_like, float
-            Background source galaxy redshift(s)
+        z_src : array_like, float, function
+            Information on the background source galaxy redshift(s). Value required depends on
+            `z_src_info` (see below).
+        z_src_info : str, optional
+            Type of redshift information provided by the `z_src` argument.
+            The following supported options are:
+
+                * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+                  It can be individual redshifts for each source galaxy when `z_src` is an
+                  array or all sources are at the same redshift when `z_src` is a float.
+
+                * 'distribution' : A redshift distribution function is provided by `z_src`.
+                  `z_src` must be a one dimentional function.
+
+                * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+                  `z_src` must be a tuple containing
+                  ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+                  the lensing efficiency and square of the lensing efficiency averaged over
+                  the galaxy redshift distribution repectively.
+
+                    .. math::
+                        \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                        {D_{L,\infty}}\right\rangle
+
+                    .. math::
+                        \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                        {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+        approx : str, optional
+            Type of computation to be made for magnifications, options are:
+
+                * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
+                  If `z_src_info='discrete'`, full computation is made for each
+                  `r_proj, z_src` pair individually. If `z_src_info='distribution'`, magnification
+                  at each value of `r_proj` is calculated as
+
+                  .. math::
+                      \mu
+                      =\left<\frac{1}{\left(1-\beta_s\kappa_{\infty}\right)^2
+                      -\left(\beta_s\gamma_{\infty}\right)^2}\right>
+                      =\frac{\int_{z_{min}}^{z_{max}}\frac{N(z)\text{d}z}
+                      {\left(1-\beta_s(z)\kappa_{\infty}\right)^2
+                      -\left(\beta_s(z)\gamma_{\infty}\right)^2}}
+                      {\int_{z_{min}}^{z_{max}} N(z)\text{d}z}
+
+                * 'order1' : Uses the weak lensing approximation of the magnification with up to
+                  first-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+                  `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+                  compute :math:`\langle \beta_s \rangle`)
+
+                  .. math::
+                      \mu \approx 1 + 2 \left<\beta_s\right>\kappa_{\infty}
+
+                * 'order2' : Uses the weak lensing approximation of the magnification with up to
+                  second-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+                  `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+                  compute :math:`\langle \beta_s \rangle`)
+
+                  .. math::
+                      \mu \approx 1 + 2 \left<\beta_s\right>\kappa_{\infty}
+                      + 3 \left<\beta_s^2\right>\kappa_{\infty}^2
+                      + \left<\beta_s^2\right>\gamma_{\infty}^2
+
+        beta_kwargs: None, dict
+            Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+            Only used if `z_src_info='distribution'`. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        verbose : bool, optional
+            If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
+            CCL backends.
 
         Returns
         -------
         mu : array_like, float
             magnification, mu.
 
-        Notes
-        -----
-        The magnification is computed taking into account just the tangential
-        shear. This is valid for spherically averaged profiles, e.g., NFW and
-        Einasto (by construction the cross shear is zero).
         """
         if self.validate_input:
             validate_argument(locals(), 'r_proj', 'float_array', argmin=0)
             validate_argument(locals(), 'z_cl', float, argmin=0)
-            validate_argument(locals(), 'z_src', 'float_array', argmin=0)
+            validate_argument(locals(), 'z_src_info', str)
+            validate_argument(locals(), 'approx', str, none_ok=True)
+            self._validate_z_src(locals())
 
         if self.halo_profile_model=='einasto' and verbose:
             print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
 
-        return self._eval_magnification(r_proj=r_proj, z_cl=z_cl, z_src=z_src)
+        if approx is None:
+            if z_src_info=='distribution':
+                core = lambda gammat, kappa: 1/((1-kappa)**2-gammat**2)
+                mu = self._pdz_weighted_avg(core, z_src, r_proj, z_cl,
+                                            integ_kwargs=beta_kwargs)
+            elif z_src_info=='discrete':
+                mu = self._eval_magnification(r_proj=r_proj, z_cl=z_cl, z_src=z_src)
+            else:
+                raise ValueError(
+                    "approx=None requires z_src_info='discrete' or 'distribution',"
+                    f"z_src_info='{z_src_info}' was provided.")
+
+        elif approx in ('order1', 'order2'):
+            beta_s_mean = self._get_beta_s_mean(
+                z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
 
-    def eval_magnification_bias(self, r_proj, z_cl, z_src, alpha):
+            kappa_inf = self._eval_convergence(r_proj, z_cl, z_src=self.z_inf)
+            gammat_inf = self._eval_tangential_shear(r_proj, z_cl, z_src=self.z_inf)
+
+            mu = 1 + 2*beta_s_mean*kappa_inf
+
+            if approx == 'order2':
+                beta_s_square_mean = self._get_beta_s_square_mean(
+                    z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
+                # Taylor expansion with up to second-order terms
+                mu += 3*beta_s_square_mean*kappa_inf**2 + beta_s_square_mean*gammat_inf**2
+
+        else:
+            raise ValueError(f"Unsupported approx (='{approx}')")
+        return mu
+
+    def eval_magnification_bias(self, r_proj, z_cl, z_src, alpha, z_src_info='discrete',
+                                approx=None, beta_kwargs=None, verbose=False):
         r"""Computes the magnification bias
 
         .. math::
@@ -771,28 +1224,143 @@ def eval_magnification_bias(self, r_proj, z_cl, z_src, alpha):
             The projected radial positions in :math:`M\!pc`.
         z_cl : float
             Galaxy cluster redshift
-        z_src : array_like, float
-            Background source galaxy redshift(s)
+        z_src : array_like, float, function
+            Information on the background source galaxy redshift(s). Value required depends on
+            `z_src_info` (see below).
         alpha : float
             Slope of the cummulative number count of background sources at a given magnitude
+        z_src_info : str, optional
+            Type of redshift information provided by the `z_src` argument.
+            The following supported options are:
+
+                * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+                  It can be individual redshifts for each source galaxy when `z_src` is an
+                  array or all sources are at the same redshift when `z_src` is a float.
+
+                * 'distribution' : A redshift distribution function is provided by `z_src`.
+                  `z_src` must be a one dimentional function.
+
+                * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+                  `z_src` must be a tuple containing
+                  ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
+                  the lensing efficiency and square of the lensing efficiency averaged over
+                  the galaxy redshift distribution repectively.
+
+                    .. math::
+                        \langle \beta_s \rangle = \left\langle \frac{D_{LS}}{D_S}\frac{D_\infty}
+                        {D_{L,\infty}}\right\rangle
+
+                    .. math::
+                        \langle \beta_s^2 \rangle = \left\langle \left(\frac{D_{LS}}
+                        {D_S}\frac{D_\infty}{D_{L,\infty}}\right)^2 \right\rangle
+
+        approx : str, optional
+            Type of computation to be made for magnification biases, options are:
+
+                * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
+                  If `z_src_info='discrete'`, full computation is made for each
+                  `r_proj, z_src` pair individually. If `z_src_info='distribution'`, magnification
+                  bias at each value of `r_proj` is calculated as
+
+                  .. math::
+                      \mu^{\alpha-1}
+                      &=\left(\left<\frac{1}{\left(1-\beta_s\kappa_{\infty}\right)^2
+                      -\left(\beta_s\gamma_{\infty}\right)^2}\right>\right)^{\alpha-1}
+                      \\\\
+                      &=\frac{\int_{z_{min}}^{z_{max}}\frac{N(z)\text{d}z}
+                      {\left(\left(1-\beta_s(z)\kappa_{\infty}\right)^2
+                      -\left(\beta_s(z)\gamma_{\infty}\right)^2\right)^{\alpha-1}}}
+                      {\int_{z_{min}}^{z_{max}} N(z)\text{d}z}
+
+                * 'order1' : Uses the weak lensing approximation of the magnification bias with up
+                  to first-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+                  `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+                  compute :math:`\langle \beta_s \rangle`)
+
+                  .. math::
+                      \mu^{\alpha-1} \approx
+                      1 + \left(\alpha-1\right)\left(2 \left<\beta_s\right>\kappa_{\infty}\right)
+
+                * 'order2' : Uses the weak lensing approximation of the magnification bias with up
+                  to second-order terms in :math:`\kappa_{\infty}` or :math:`\gamma_{\infty}`.
+                  `z_src_info` must be either 'beta', or 'distribution' (that will be used to
+                  compute :math:`\langle \beta_s \rangle`)
+
+                  .. math::
+                      \mu^{\alpha-1} \approx
+                      1 &+ \left(\alpha-1\right)\left(2 \left<\beta_s\right>\kappa_{\infty}\right)
+                      \\\\
+                      &+ \left(\alpha-1\right)\left(\left<\beta_s^2\right>\gamma_{\infty}^2\right)
+                      \\\\
+                      &+ \left(2\alpha-1\right)\left(\alpha-1\right)
+                      \left(\left<\beta_s^2\right>\kappa_{\infty}^2\right)
+
+        beta_kwargs: None, dict
+            Extra arguments for the `compute_beta_s_mean, compute_beta_s_square_mean` functions.
+            Only used if `z_src_info='distribution'`. Possible keys are:
+
+                * 'zmin' (None, float) : Minimum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=None)
+                * 'zmax' (float) : Maximum redshift to be set as the source of the galaxy
+                  when performing the sum. (default=10.0)
+                * 'delta_z_cut' (float) : Redshift cut so that `zmin` = `z_cl` + `delta_z_cut`.
+                  `delta_z_cut` is ignored if `z_min` is already provided. (default=0.1)
+
+        verbose : bool, optional
+            If True, the Einasto slope (alpha_ein) is printed out. Only availble for the NC and
+            CCL backends.
 
         Returns
         -------
         mu_bias : array_like, float
             magnification bias.
 
-        Notes
-        -----
-        The magnification is computed taking into account just the tangential
-        shear. This is valid for spherically averaged profiles, e.g., NFW and
-        Einasto (by construction the cross shear is zero).
         """
         if self.validate_input:
             validate_argument(locals(), 'r_proj', 'float_array', argmin=0)
             validate_argument(locals(), 'z_cl', float, argmin=0)
-            validate_argument(locals(), 'z_src', 'float_array', argmin=0)
+            validate_argument(locals(), 'z_src_info', str)
             validate_argument(locals(), 'alpha', 'float_array')
-        return self._eval_magnification_bias(r_proj=r_proj, z_cl=z_cl, z_src=z_src, alpha=alpha)
+            validate_argument(locals(), 'approx', str, none_ok=True)
+            self._validate_z_src(locals())
+
+        if self.halo_profile_model=='einasto' and verbose:
+            print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
+
+        if approx is None:
+            # z_src (float or array) is redshift
+            if z_src_info=='distribution':
+                core = lambda gammat, kappa: 1/((1-kappa)**2-gammat**2)**(alpha-1)
+                mu_bias = self._pdz_weighted_avg(core, z_src, r_proj, z_cl,
+                                                 integ_kwargs=beta_kwargs)
+            elif z_src_info=='discrete':
+                mu_bias = self._eval_magnification_bias(
+                    r_proj=r_proj, z_cl=z_cl, z_src=z_src, alpha=alpha)
+            else:
+                raise ValueError(
+                    "approx=None requires z_src_info='discrete' or 'distribution',"
+                    f"z_src_info='{z_src_info}' was provided.")
+
+        elif approx in ('order1', 'order2'):
+            beta_s_mean = self._get_beta_s_mean(
+                z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
+
+            kappa_inf = self._eval_convergence(r_proj, z_cl, z_src=self.z_inf)
+            gammat_inf = self._eval_tangential_shear(r_proj, z_cl, z_src=self.z_inf)
+
+            mu_bias = 1 + (alpha-1)*(2*beta_s_mean*kappa_inf)
+
+            if approx == 'order2':
+                beta_s_square_mean = self._get_beta_s_square_mean(
+                    z_cl, z_src, z_src_info=z_src_info, beta_kwargs=beta_kwargs)
+                # Taylor expansion with up to second-order terms
+                mu_bias += (alpha-1)*(beta_s_square_mean*gammat_inf**2)\
+                           +(2*alpha-1)*(alpha-1)*beta_s_square_mean*kappa_inf**2
+
+        else:
+            raise ValueError(f"Unsupported approx (='{approx}')")
+
+        return mu_bias
 
     def eval_rdelta(self, z_cl):
         r"""Retrieves the radius for mdelta
@@ -905,3 +1473,39 @@ def convert_mass_concentration(self, z_cl, massdef=None, delta_mdef=None,
 
         return self._convert_mass_concentration(z_cl, massdef, delta_mdef,
                                                 halo_profile_model, alpha)
+
+    def _validate_z_src(self, loc_dict):
+        r"""Validation for z_src according to z_src_info. The conditions are:
+
+            * z_src_info='discrete' : z_src must be array or float.
+            * z_src_info='distribution' : z_src must be a one dimentional function.
+            * z_src_info='beta' : z_src must be a tuple containing
+              ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`).
+
+        Also, if approx is provided and not None, z_src_info must be 'distribution' or 'beta'.
+
+        Parameters
+        ----------
+        locals_dict: dict
+            Should be the call locals()
+        """
+        if loc_dict['z_src_info']=='discrete':
+            validate_argument(loc_dict, 'z_src', 'float_array', argmin=0)
+        elif loc_dict['z_src_info']=='distribution':
+            validate_argument(loc_dict, 'z_src', 'function', none_ok=False)
+            beta_kwargs = {} if loc_dict['beta_kwargs'] is None else loc_dict['beta_kwargs']
+            _def_keys = ['zmin', 'zmax', 'delta_z_cut']
+            if any(key not in _def_keys for key in beta_kwargs):
+                raise KeyError(f'beta_kwargs must contain only {_def_keys} keys, '
+                               f' {beta_kwargs.keys()} provided.')
+        elif loc_dict['z_src_info']=='beta':
+            validate_argument(loc_dict, 'z_src', 'array')
+            beta_info = {'beta_s_mean':loc_dict['z_src'][0],
+                         'beta_s_square_mean':loc_dict['z_src'][1]}
+            validate_argument(beta_info, 'beta_s_mean', 'float_array')
+            validate_argument(beta_info, 'beta_s_square_mean', 'float_array')
+        if loc_dict.get('approx') and loc_dict['z_src_info'] not in ('distribution', 'beta'):
+            approx, z_src_info = loc_dict['approx'], loc_dict['z_src_info']
+            raise ValueError(
+                f"approx='{approx}' requires z_src_info='distribution' or 'beta', "
+                f"z_src_info='{z_src_info}' was provided.")
diff --git a/clmm/utils.py b/clmm/utils.py
index 9c2e62f66..4621c733b 100644
--- a/clmm/utils.py
+++ b/clmm/utils.py
@@ -3,10 +3,10 @@
 import numpy as np
 from astropy import units as u
 from scipy.stats import binned_statistic
-from scipy.special import gamma, gammainc
 from scipy.integrate import quad, cumulative_trapezoid, simps
 from scipy.interpolate import interp1d
 from .constants import Constants as const
+from . import z_distributions as zdist
 
 
 def compute_nfw_boost(rvals, rs=1000, b0=0.1) :
@@ -487,7 +487,8 @@ def _patch_rho_crit_to_cd2018(rho_crit_external):
     float: (float, int, np.floating, np.integer),
     int: (int, np.integer),
     'float_array': (float, int, np.floating, np.integer),
-    'int_array': (int, np.integer)
+    'int_array': (int, np.integer),
+    'array': (list, tuple, np.ndarray),
     }
 
 def _is_valid(arg, valid_type):
@@ -508,6 +509,8 @@ def _is_valid(arg, valid_type):
     valid: bool
         Is argument valid
     """
+    if valid_type=='function':
+        return callable(arg)
     return (isinstance(arg[0], _valid_types[valid_type])
                 if (valid_type in ('int_array', 'float_array') and np.iterable(arg))
                 else isinstance(arg, _valid_types.get(valid_type, valid_type)))
@@ -572,7 +575,7 @@ def validate_argument(loc, argname, valid_type, none_ok=False, argmin=None, argm
 
 def _integ_pzfuncs(pzpdf, pzbins, zmin=0., zmax=5, kernel=lambda z: 1., is_unique_pzbins=False, ngrid=1000):
     r"""
-    Integrates the product of a photo-z pdf with a given kernel.
+    Integrates the product of a photo-z pdf with a given kernel. 
     This function was created to allow for data with different photo-z binnings.
 
     Parameters
@@ -651,7 +654,7 @@ def compute_for_good_redshifts(function, z1, z2, bad_value, error_message):
         res = function(z1, z2)
     return res
 
-def compute_beta(z_s, z_cl, cosmo):
+def compute_beta(z_src, z_cl, cosmo):
     r"""Geometric lensing efficicency
 
     .. math::
@@ -661,10 +664,10 @@ def compute_beta(z_s, z_cl, cosmo):
 
     Parameters
     ----------
+    z_src:  float
+        Source galaxy redshift
     z_cl: float
-            Galaxy cluster redshift
-    z_s:  float
-            Source galaxy  redshift
+        Galaxy cluster redshift
     cosmo: clmm.Cosmology
         CLMM Cosmology object
 
@@ -673,23 +676,23 @@ def compute_beta(z_s, z_cl, cosmo):
     float
         Geometric lensing efficicency
     """
-    beta = np.heaviside(z_s-z_cl, 0) * cosmo.eval_da_z1z2(z_cl, z_s) / cosmo.eval_da(z_s)
+    beta = np.heaviside(z_src-z_cl, 0) * cosmo.eval_da_z1z2(z_cl, z_src) / cosmo.eval_da(z_src)
     return beta
 
-def compute_beta_s(z_s, z_cl, z_inf, cosmo):
+def compute_beta_s(z_src, z_cl, z_inf, cosmo):
     r"""Geometric lensing efficicency ratio
 
     .. math::
-        \beta_s = \beta(z_s)/\beta(z_{inf})
+        \beta_s = \beta(z_{src})/\beta(z_{inf})
 
     Parameters
     ----------
+    z_src: float
+        Source galaxy redshift
     z_cl: float
-            Galaxy cluster redshift
-    z_s: float
-            Source galaxy redshift
+        Galaxy cluster redshift
     z_inf: float
-            Redshift at infinity
+        Redshift at infinity
     cosmo: clmm.Cosmology
         CLMM Cosmology object
 
@@ -698,15 +701,48 @@ def compute_beta_s(z_s, z_cl, z_inf, cosmo):
     float
         Geometric lensing efficicency ratio
     """
-    beta_s = compute_beta(z_s, z_cl, cosmo) / compute_beta(z_inf, z_cl, cosmo)
+    beta_s = compute_beta(z_src, z_cl, cosmo) / compute_beta(z_inf, z_cl, cosmo)
     return beta_s
 
+def compute_beta_s_func(z_src, z_cl, z_inf, cosmo, func, *args, **kwargs):
+    r"""Geometric lensing efficicency ratio times a value of a function
+
+    .. math::
+        \beta_{s}\times \text{func} = \beta_s(z_{src}, z_{cl}, z_{inf})
+        \times\text{func}(*args,\ **kwargs)
+
+    Parameters
+    ----------
+    z_src: float
+        Source galaxy redshift
+    z_cl: float
+        Galaxy cluster redshift
+    z_inf: float
+        Redshift at infinity
+    cosmo: clmm.Cosmology
+        CLMM Cosmology object
+    func: callable
+        A scalar function
+    *args: positional arguments
+        args to be passed to `func`
+    **kwargs: keyword arguments
+        kwargs to be passed to `func`
+
+    Returns
+    -------
+    float
+        Geometric lensing efficicency ratio
+    """
+    beta_s = compute_beta(z_src, z_cl, cosmo) / compute_beta(z_inf, z_cl, cosmo)
+    beta_s_func = beta_s * func(*args, **kwargs)
+    return beta_s_func
+
 def compute_beta_mean(z_cl, cosmo, zmax=10.0, delta_z_cut=0.1, zmin=None, z_distrib_func=None):
     r"""Mean value of the geometric lensing efficicency
 
     .. math::
-       \left<\beta\right> = \frac{\sum_{z = z_{min}}^{z_{max}}\beta(z)p(z)}
-       {\sum_{z = z_{min}}^{z_{max}}p(z)}
+       \left<\beta\right> = \frac{\int_{z = z_{min}}^{z_{max}}\beta(z)N(z)}
+       {\int_{z = z_{min}}^{z_{max}}N(z)}
 
     Parameters
     ----------
@@ -731,8 +767,8 @@ def compute_beta_mean(z_cl, cosmo, zmax=10.0, delta_z_cut=0.1, zmin=None, z_dist
     float
         Mean value of the geometric lensing efficicency
     """
-    if z_distrib_func == None:
-        z_distrib_func = _chang_z_distrib
+    if z_distrib_func is None:
+        z_distrib_func = zdist.chang2013
     def integrand(z_i, z_cl=z_cl, cosmo=cosmo):
         return compute_beta(z_i, z_cl, cosmo) * z_distrib_func(z_i)
 
@@ -746,8 +782,8 @@ def compute_beta_s_mean(z_cl, z_inf, cosmo, zmax=10.0, delta_z_cut=0.1, zmin=Non
     r"""Mean value of the geometric lensing efficicency ratio
 
     .. math::
-       \left<\beta_s\right> =\frac{\sum_{z = z_{min}}^{z_{max}}\beta_s(z)p(z)}
-       {\sum_{z = z_{min}}^{z_{max}}p(z)}
+       \left<\beta_s\right> =\frac{\int_{z = z_{min}}^{z_{max}}\beta_s(z)N(z)}
+       {\int_{z = z_{min}}^{z_{max}}N(z)}
 
     Parameters
     ----------
@@ -763,22 +799,19 @@ def compute_beta_s_mean(z_cl, z_inf, cosmo, zmax=10.0, delta_z_cut=0.1, zmin=Non
     delta_z_cut: float, optional
         Redshift interval to be summed with :math:`z_{cl}` to return :math:`z_{min}`.
         This feature is not used if :math:`z_{min}` is provided by the user. Default: 0.1
-    z_distrib_func: one-parameter function
-        Redshift distribution function. Default is Chang et al (2013) distribution function.
     zmin: float, None, optional
         Minimum redshift to be set as the source of the galaxy when performing the sum.
         Default: None
     z_distrib_func: one-parameter function, optional
         Redshift distribution function. Default is Chang et al (2013) distribution function.
 
-
     Returns
     -------
     float
         Mean value of the geometric lensing efficicency ratio
     """
-    if z_distrib_func == None:
-        z_distrib_func = _chang_z_distrib
+    if z_distrib_func is None:
+        z_distrib_func = zdist.chang2013
 
     def integrand(z_i, z_cl=z_cl, z_inf=z_inf, cosmo=cosmo):
         return compute_beta_s(z_i, z_cl, z_inf, cosmo) * z_distrib_func(z_i)
@@ -792,15 +825,15 @@ def compute_beta_s_square_mean(z_cl, z_inf, cosmo, zmax=10.0, delta_z_cut=0.1, z
     r"""Mean square value of the geometric lensing efficiency ratio
 
     .. math::
-       \left<\beta_s^2\right> =\frac{\sum_{z = z_{min}}^{z_{max}}\beta_s^2(z)p(z)}
-       {\sum_{z = z_{min}}^{z_{max}}p(z)}
+       \left<\beta_s^2\right> =\frac{\int_{z = z_{min}}^{z_{max}}\beta_s^2(z)N(z)}
+       {\int_{z = z_{min}}^{z_{max}}N(z)}
 
     Parameters
     ----------
     z_cl: float
-            Galaxy cluster redshift
+        Galaxy cluster redshift
     z_inf: float
-            Redshift at infinity
+        Redshift at infinity
     cosmo: clmm.Cosmology
         CLMM Cosmology object
     zmax: float
@@ -820,8 +853,8 @@ def compute_beta_s_square_mean(z_cl, z_inf, cosmo, zmax=10.0, delta_z_cut=0.1, z
     float
         Mean square value of the geometric lensing efficicency ratio.
     """
-    if z_distrib_func == None:
-        z_distrib_func = _chang_z_distrib
+    if z_distrib_func is None:
+        z_distrib_func = zdist.chang2013
 
     def integrand(z_i, z_cl=z_cl, z_inf=z_inf, cosmo=cosmo):
         return compute_beta_s(z_i, z_cl, z_inf, cosmo)**2 * z_distrib_func(z_i)
@@ -831,49 +864,6 @@ def integrand(z_i, z_cl=z_cl, z_inf=z_inf, cosmo=cosmo):
     Bs_square_mean = quad(integrand, zmin, zmax)[0] / quad(z_distrib_func, zmin, zmax)[0]
     return Bs_square_mean
 
-def _chang_z_distrib(redshift, is_cdf=False):
-    """
-    A private function that returns the Chang et al (2013) unnormalized galaxy redshift
-    distribution function, with the fiducial set of parameters.
-
-    Parameters
-    ----------
-    redshift : float
-        Galaxy redshift
-    is_cdf : bool, optional
-        If True, returns cumulative distribution function. Default: False
-
-    Returns
-    -------
-    The value of the distribution at z
-    """
-    alpha, beta, redshift0 = 1.24, 1.01, 0.51
-    if is_cdf:
-        return redshift0**(alpha+1)*gammainc((alpha+1)/beta, (redshift/redshift0)**beta)/beta*gamma((alpha+1)/beta)
-    else:
-        return (redshift**alpha)*np.exp(-(redshift/redshift0)**beta)
-
-def _srd_z_distrib(redshift, is_cdf=False):
-    """
-    A private function that returns the unnormalized galaxy redshift distribution function used in
-    the LSST/DESC Science Requirement Document (arxiv:1809.01669).
-
-    Parameters
-    ----------
-    redshift : float
-        Galaxy redshift
-    is_cdf : bool, optional
-        If True, returns cumulative distribution function. Default: False
-
-    Returns
-    -------
-    The value of the distribution at z
-    """
-    alpha, beta, redshift0 = 2., 0.9, 0.28
-    if is_cdf:
-        return redshift0**(alpha+1)*gammainc((alpha+1)/beta, (redshift/redshift0)**beta)/beta*gamma((alpha+1)/beta)
-    else:
-        return (redshift**alpha)*np.exp(-(redshift/redshift0)**beta)
 
 def _draw_random_points_from_distribution(xmin, xmax, nobj, dist_func, xstep=0.001):
     """Draw random points with a given distribution.
diff --git a/clmm/z_distributions.py b/clmm/z_distributions.py
new file mode 100644
index 000000000..3754280d8
--- /dev/null
+++ b/clmm/z_distributions.py
@@ -0,0 +1,70 @@
+"""Redshift distribution functions"""
+import numpy as np
+from scipy.special import gamma, gammainc
+
+def _functional_form(redshift, alpha, beta, redshift0, is_cdf=False):
+    """
+    A private function that returns the functionnal form of the redshift distribution used in
+    Chang et al (2013):
+
+    .. math::
+       P(z) = z^{\\alpha}\times\exp^{\left(-\frac{z}{z0}^\beta\right)}
+
+    Parameters
+    ----------
+    redshift : float
+        Galaxy redshift
+    alpha, beta, z0 : floats
+        Parameters describing the function
+    is_cdf : bool
+        If True, returns cumulative function.
+
+    Returns
+    -------
+    The value of the function at z
+    """
+    if is_cdf:
+        return redshift0**(alpha+1)*gammainc((alpha+1)/beta, (redshift/redshift0)**beta)\
+    /beta*gamma((alpha+1)/beta)
+    else:
+        return (redshift**alpha)*np.exp(-(redshift/redshift0)**beta)
+
+
+def chang2013(redshift, is_cdf=False):
+    """
+    Chang et al (2013) unnormalized galaxy redshift distribution function, with the fiducial
+    set of parameters.
+
+    Parameters
+    ----------
+    redshift : float
+        Galaxy redshift
+    is_cdf : bool
+        If True, returns cumulative distribution function.
+
+    Returns
+    -------
+    The value of the distribution at z
+    """
+    alpha, beta, redshift0 = 1.24, 1.01, 0.51
+    return _functional_form(redshift, alpha, beta, redshift0, is_cdf)
+
+
+def desc_srd(redshift, is_cdf=False):
+    """
+    Unnormalized galaxy redshift distribution function used in
+    the LSST/DESC Science Requirement Document (arxiv:1809.01669).
+
+    Parameters
+    ----------
+    redshift : float
+        Galaxy redshift
+    is_cdf : bool
+        If True, returns cumulative distribution function.
+
+    Returns
+    -------
+    The value of the distribution at z
+    """
+    alpha, beta, redshift0 = 2., 0.9, 0.28
+    return _functional_form(redshift, alpha, beta, redshift0, is_cdf)
diff --git a/docs/doc-config.ini b/docs/doc-config.ini
index 495220f48..dd3f60c29 100644
--- a/docs/doc-config.ini
+++ b/docs/doc-config.ini
@@ -8,6 +8,7 @@ gcdata
 theory
 plotting
 utils
+z_distributions
 support
 
 DEMO
diff --git a/examples/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb b/examples/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb
index a018b2aff..3f1278327 100644
--- a/examples/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb
+++ b/examples/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb
@@ -12,7 +12,7 @@
     "This notebook demonstrates how to use `clmm` to estimate a WL halo mass from observations of a galaxy cluster when source galaxies follow a given distribution (the Chang et al. (2013) implemented in `clmm`). It uses several functionalities of the support `mock_data` module to produce mock datasets.\n",
     "\n",
     "- Setting things up, with the proper imports.\n",
-    "- Generating 2 datasets: an ideal datasetwith source galaxies following the Chang et al. (2013) redshift distribution; a noisy dataset where photoz errors and shape noise are also included. \n",
+    "- Generating 2 datasets: an ideal dataset with source galaxies following the Chang et al. (2013) redshift distribution; a noisy dataset where photoz errors and shape noise are also included. \n",
     "- Computing the binned reduced tangential shear profile, for the 2 datasets, using logarithmic binning.\n",
     "- Setting up the \"single source plane\" model (model1) and a model accounting for the redshift distribution (model2). As already seen in Example2, model1 will yield a biased mass reconstruction. Accounting for the redshift distribution in the model (model2) solves that issue. \n",
     "- We also add two models using the methods (to use the statistics of angular diameter distance ratio) from Applegate et al. (2014) and from Schrabback et al. (2018) respectively.\n",
@@ -388,8 +388,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Note \"data\" can be skipped here.\n",
-    "def model_reduced_tangential_shear_zdistrib(radius, logm, data, catalog, profile): \n",
+    "def model_reduced_tangential_shear_zdistrib(radius, logm, catalog, profile): \n",
     "    m = 10**logm\n",
     "    gt_model = []\n",
     "    for i in range(len(radius)):\n",
@@ -409,7 +408,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "- Third, we consider two models that use the statistics of angular diameter distance ratio ($\\beta_s$, i.e. $<\\beta_s>$, $<\\beta^2_s>$). They are similar to the ones we showed in the notebook for paper v1.0 (`gt_and_use_case.ipynb`) and the notebook for real data (Example4)."
+    "- Third, we consider two models that use the statistics of angular diameter distance ratio ($\\beta_s$, i.e. $\\langle\\beta_s\\rangle$, $\\langle\\beta^2_s\\rangle>$). They are similar to the ones we showed in the notebook for paper v1.0 (`gt_and_use_case.ipynb`) and the notebook for real data (Example4)."
    ]
   },
   {
@@ -427,12 +426,10 @@
     "# clmm/theory/func_layer.py: compute_reduced_tangential_shear() \n",
     "# clmm/theory/parent_class.py: eval_reduced_tangential_shear()\n",
     "\n",
-    "z_inf = 1000.\n",
-    "\n",
-    "\n",
-    "# The method from Applegate et al. (2014).\n",
-    "# The argument \"profile\" is for r; if r is defined outside, we can use r as an argument instead.\n",
+    "# The method from Applegate et al. (2014), called 'order1' in CLMM.\n",
+    "# The argument \"profile\" is only used to get the radius values r = profile['radius']\n",
     "# The default mass is M200m; use massdef='critical' for M200c.\n",
+    "\n",
     "def model_reduced_tangential_shear_applegate14(logm, catalog, profile):\n",
     "    z_values = catalog.galcat['z']\n",
     "    z_cl = cluster_z\n",
@@ -441,19 +438,19 @@
     "            mdelta=10**logm, # Mass of the cluster [M_sun]\n",
     "            cdelta=concentration, # Concentration of the cluster\n",
     "            z_cluster=cluster_z, # Redshift of the cluster\n",
-    "            z_source=z_inf, # Redshift value can be arbitrary because it's not used here\n",
+    "            z_source=clmm.z_distributions.chang2013, # We provide the redshift distribution (default: Chang et al. 2013) for calculating the beta_s statistics\n",
     "            cosmo=cosmo,\n",
     "            delta_mdef=200,\n",
     "            #massdef='critical',\n",
     "            halo_profile_model='nfw',\n",
-    "            z_src_model='applegate14', \n",
+    "            z_src_info='distribution',\n",
+    "            approx='order1',\n",
     "            #beta_s_mean=None, # beta_s_mean and beta_s_square_mean can be provided (default: None)\n",
-    "            #beta_s_square_mean=None,\n",
-    "            z_distrib_func=clmm.utils._chang_z_distrib, # We provide the redshift distribution (default: Chang et al. 2013) for calculating the beta_s statistics  \n",
+    "            #beta_s_square_mean=None \n",
     "    )\n",
     "    return gt_model\n",
     "\n",
-    "# Similarly, we also consider the method from Schrabback et al. (2018).\n",
+    "# Similarly, we also consider the method from Schrabback et al. (2018), called 'order12' in CLMM.\n",
     "def model_reduced_tangential_shear_schrabback18(logm, catalog, profile):\n",
     "    z_values = catalog.galcat['z']\n",
     "    gt_model = clmm.compute_reduced_tangential_shear(\n",
@@ -461,15 +458,15 @@
     "            mdelta=10**logm, # Mass of the cluster [M_sun]\n",
     "            cdelta=concentration, # Concentration of the cluster\n",
     "            z_cluster=cluster_z, # Redshift of the cluster\n",
-    "            z_source=z_inf, # Redshift value can be arbitrary because it's not used here\n",
+    "            z_source=clmm.z_distributions.chang2013, # We provide the redshift distribution (default: Chang et al. 2013) for calculating the beta_s statistics\n",
     "            cosmo=cosmo,\n",
     "            delta_mdef=200,\n",
     "            #massdef='critical',\n",
     "            halo_profile_model='nfw',\n",
-    "            z_src_model='schrabback18',\n",
+    "            z_src_info='distribution',\n",
+    "            approx='order2',\n",
     "            #beta_s_mean=None, \n",
-    "            #beta_s_square_mean=None, \n",
-    "            z_distrib_func=clmm.utils._chang_z_distrib,\n",
+    "            #beta_s_square_mean=None\n",
     "    )\n",
     "    return gt_model\n",
     "\n"
@@ -492,13 +489,13 @@
     "r = cl_ideal.profile['radius']\n",
     "\n",
     "gt_model_ideal_singlez = model_reduced_tangential_shear_singlez(r, logm, cl_ideal.profile['z'])\n",
-    "gt_model_ideal_zdistrib = model_reduced_tangential_shear_zdistrib(r, logm, ideal_data, cl_ideal, cl_ideal.profile)\n",
+    "gt_model_ideal_zdistrib = model_reduced_tangential_shear_zdistrib(r, logm, cl_ideal, cl_ideal.profile)\n",
     "gt_model_ideal_applegate14 = model_reduced_tangential_shear_applegate14(logm, cl_ideal, cl_ideal.profile)\n",
     "gt_model_ideal_schrabback18 = model_reduced_tangential_shear_schrabback18(logm, cl_ideal, cl_ideal.profile)\n",
     "\n",
     "# Noisy data can have some source redshifts that are lower than the cluster redshift.\n",
     "gt_model_noisy_singlez = model_reduced_tangential_shear_singlez(r,logm, cl_noisy.profile['z'])\n",
-    "gt_model_noisy_zdistrib = model_reduced_tangential_shear_zdistrib(r,logm, noisy_data, cl_noisy, cl_noisy.profile)\n",
+    "gt_model_noisy_zdistrib = model_reduced_tangential_shear_zdistrib(r,logm, cl_noisy, cl_noisy.profile)\n",
     "gt_model_noisy_applegate14 = model_reduced_tangential_shear_applegate14(logm, cl_noisy, cl_noisy.profile)\n",
     "gt_model_noisy_schrabback18 = model_reduced_tangential_shear_schrabback18(logm, cl_noisy, cl_noisy.profile)"
    ]
@@ -581,7 +578,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "popt,pcov = fitters['curve_fit'](lambda r, logm:model_reduced_tangential_shear_zdistrib(r, logm, ideal_data, cl_ideal, cl_ideal.profile), \n",
+    "popt,pcov = fitters['curve_fit'](lambda r, logm:model_reduced_tangential_shear_zdistrib(r, logm, cl_ideal, cl_ideal.profile), \n",
     "                        cl_ideal.profile['radius'], \n",
     "                        cl_ideal.profile['gt'], \n",
     "                        cl_ideal.profile['gt_err'], bounds=[10.,16.])\n",
@@ -620,7 +617,7 @@
     "\n",
     "\n",
     "\n",
-    "popt,pcov = fitters['curve_fit'](lambda r, logm:model_reduced_tangential_shear_zdistrib(r, logm, noisy_data, cl_noisy, cl_noisy.profile), \n",
+    "popt,pcov = fitters['curve_fit'](lambda r, logm:model_reduced_tangential_shear_zdistrib(r, logm, cl_noisy, cl_noisy.profile), \n",
     "                        cl_noisy.profile['radius'], \n",
     "                        cl_noisy.profile['gt'], \n",
     "                        cl_noisy.profile['gt_err'], bounds=[10.,16.])\n",
@@ -710,8 +707,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gt_est_ideal_zdistrib = model_reduced_tangential_shear_zdistrib(r,np.log(m_est_ideal_zdistrib)/np.log(10), ideal_data, cl_ideal, cl_ideal.profile)\n",
-    "gt_est_noisy_zdistrib = model_reduced_tangential_shear_zdistrib(r,np.log(m_est_noisy_zdistrib)/np.log(10), noisy_data, cl_noisy, cl_noisy.profile)\n",
+    "gt_est_ideal_zdistrib = model_reduced_tangential_shear_zdistrib(r,np.log(m_est_ideal_zdistrib)/np.log(10), cl_ideal, cl_ideal.profile)\n",
+    "gt_est_noisy_zdistrib = model_reduced_tangential_shear_zdistrib(r,np.log(m_est_noisy_zdistrib)/np.log(10), cl_noisy, cl_noisy.profile)\n",
     "gt_est_ideal_applegate14 = model_reduced_tangential_shear_applegate14(np.log(m_est_ideal_applegate14)/np.log(10), cl_ideal, cl_ideal.profile)\n",
     "gt_est_noisy_applegate14 = model_reduced_tangential_shear_applegate14(np.log(m_est_noisy_applegate14)/np.log(10), cl_noisy, cl_noisy.profile)\n",
     "gt_est_ideal_schrabback18 = model_reduced_tangential_shear_schrabback18(np.log(m_est_ideal_schrabback18)/np.log(10), cl_ideal, cl_ideal.profile)\n",
@@ -783,9 +780,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "conda-clmmenv",
+   "display_name": "wrk",
    "language": "python",
-   "name": "conda-clmmenv"
+   "name": "wrk"
   },
   "language_info": {
    "codemirror_mode": {
@@ -797,7 +794,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.0"
+   "version": "3.10.6"
   }
  },
  "nbformat": 4,
diff --git a/examples/Example_different_redshift_input_type.ipynb b/examples/Example_different_redshift_input_type.ipynb
new file mode 100644
index 000000000..010b58c87
--- /dev/null
+++ b/examples/Example_different_redshift_input_type.ipynb
@@ -0,0 +1,730 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "# Model profiles using different type of source redshift information as input"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this example we model lensing profiles by giving as input either : \n",
+    "- discrete source redshifts, \n",
+    "- a redshift distribution function,\n",
+    "- the value of the mean beta parameters : \n",
+    "$\\langle \\beta_s \\rangle = \\left\\langle \\frac{D_{LS}}{D_S}\\frac{D_\\infty}{D_{L,\\infty}}\\right\\rangle$ ,\n",
+    "$\\langle \\beta_s^2 \\rangle = \\left\\langle \\left(\\frac{D_{LS}}{D_S}\\frac{D_\\infty}{D_{L,\\infty}}\\right)^2 \\right\\rangle$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\", message='.*(!).*')\n",
+    "import os\n",
+    "## Uncomment the following line if you want to use a specific modeling backend among 'ct' (cluster-toolkit), 'ccl' (CCL) or 'nc' (Numcosmo). Default is 'ccl'\n",
+    "#os.environ['CLMM_MODELING_BACKEND'] = 'nc'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try: import clmm\n",
+    "except:\n",
+    "    import notebook_install\n",
+    "    notebook_install.install_clmm_pipeline(upgrade=False)\n",
+    "    import clmm\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib\n",
+    "matplotlib.rcParams.update({'font.size': 16})\n",
+    "matplotlib.rcParams.update({'figure.figsize': (20,6)})"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Make sure we know which version we're using"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "clmm.__version__"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Import mock data module and setup the configuration "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from clmm.support import mock_data as mock\n",
+    "from clmm import Cosmology\n",
+    "from clmm.z_distributions import *"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Mock data generation requires a defined cosmology"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mock_cosmo = Cosmology(H0 = 70.0, Omega_dm0 = 0.27 - 0.045, Omega_b0 = 0.045, Omega_k0 = 0.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Mock data generation requires some cluster information"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cosmo = mock_cosmo\n",
+    "\n",
+    "# cluster properties from https://arxiv.org/pdf/1611.03866.pdf\n",
+    "cluster_id = \"SPT-CL J0000−5748\"\n",
+    "cluster_m = 4.56e14 # M500,c\n",
+    "cluster_z = 0.702\n",
+    "cluster_ra = 0.2499 \n",
+    "cluster_dec = -57.8064\n",
+    "concentration = 5 # (arbitrary value, not from the paper)\n",
+    "\n",
+    "#source redshift distribution properties\n",
+    "cluster_beta_s_mean = 0.466\n",
+    "cluster_beta_s2_mean = 0.243\n",
+    "ngal_density= 26. *100 # density of source galaxies per arcmin^2 # (arbitrary value, not from the paper)\n",
+    "model_z_distrib_dict = {'func':desc_srd, 'name':\"desc_srd\"}\n",
+    "delta_z_cut = 0.1\n",
+    "zsrc_min = cluster_z + delta_z_cut\n",
+    "zsrc_max = 3.0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 1 - Defining the different inputs for the source redshifts"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Discrete redshifts"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Generate the mock source catalog\n",
+    "\n",
+    "The CLMM mock data generation will provide, among other things, a redshift value for each background galaxy that is draw from the redshift distribution given by `model_z_distrib_dict`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(0)\n",
+    "source_catalog = mock.generate_galaxy_catalog(cluster_m, cluster_z, concentration, cosmo, model_z_distrib_dict['name'], delta_so=500,\n",
+    "                                              massdef='critical', zsrc_min=zsrc_min, zsrc_max=zsrc_max, ngal_density=ngal_density, \n",
+    "                                              cluster_ra=cluster_ra, cluster_dec=cluster_dec)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Beta parameters\n",
+    "From this udnerlying redshift distribution, one may directly compute the average $\\langle\\beta_s\\rangle$ and $\\langle\\beta_s^2\\rangle$ quantities"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z_inf = 1000\n",
+    "\n",
+    "beta_s_mean        = clmm.utils.compute_beta_s_mean(cluster_z, z_inf, cosmo, zmax=zsrc_max, delta_z_cut=delta_z_cut, zmin=None, z_distrib_func=model_z_distrib_dict['func'])\n",
+    "beta_s_square_mean = clmm.utils.compute_beta_s_square_mean(cluster_z, z_inf, cosmo, zmax=zsrc_max, delta_z_cut=delta_z_cut, zmin=None, z_distrib_func=model_z_distrib_dict['func'])\n",
+    "\n",
+    "print (\"$<\\\\beta_s>$ = \", round(beta_s_mean,3), \",\\n$<\\\\beta_s^2>$ = \", round(beta_s_square_mean,3))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualisation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = np.linspace(0,zsrc_max,1000)\n",
+    "\n",
+    "plt.hist(source_catalog['z'], bins=50, alpha=0.3, density=True, label='mock source redshift');\n",
+    "plt.axvline(zsrc_min, color='red', label='requested zmin')\n",
+    "plt.text(4, 0.6, '$\\\\langle\\\\beta_s\\\\rangle$ = '+str(round(beta_s_mean,3)) + ', $\\\\langle\\\\beta_s^2\\\\rangle$ = ' + str(round(beta_s_square_mean,3)))\n",
+    "#here we multiply by a constant for visualisation purposes\n",
+    "plt.plot(z,model_z_distrib_dict['func'](z)*25, linestyle='dashed',label='redshift distribution fct (arbitrary normalisation)')\n",
+    "plt.xlabel('$z_s$')\n",
+    "plt.ylabel('pdf')\n",
+    "plt.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here we use a zmax=3 (different from default values) to highlight the importance of specifying the zmax when computing the modeling. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 2 - Compute models"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here we are computing the models for the tangential shear, reduced shear, convergence, magnification and magnification bias."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Profile from mock data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First we compute the cluster profile based on the mock source catalog produced in the previous section."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gc_object = clmm.GalaxyCluster(cluster_id, cluster_ra, cluster_dec, \n",
+    "                               cluster_z, source_catalog)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gc_object.compute_tangential_and_cross_components()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gc_object.make_radial_profile('Mpc', bins=10, cosmo=cosmo)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Different ways of modeling the reduced shear"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#define radius\n",
+    "rr = np.logspace(np.log10(0.2), np.log10(5), 10)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### **Case 1 : Discrete redshift and exact formula**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### In case we know the discrete source redshift, we can compute the reduced shear for each source galaxy and take the average at a given radius. \n",
+    "#### This may take a bit of time, depending on the size of source redshift catalog and number of radius points."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "#in case we know the discrete source redshift, we can compute the reduced shear for each source galaxy and take the average at a given radius.\n",
+    "gt_discrete = np.zeros(rr.size)\n",
+    "\n",
+    "for i in range(rr.size):\n",
+    "    gt_discrete[i] = np.mean(clmm.theory.compute_reduced_tangential_shear(rr[i], cluster_m, concentration, cluster_z, \n",
+    "                                            source_catalog['z'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='discrete', approx=None))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### **Case 2 : Redshift distribution and exact formula**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### If we don't know the exact source redshift but we know the source mean redshift distribution function, we can give it as input. \\\\\n",
+    "#### In this case, we are integrating over the distribution function so it may be quite slow."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "gt_distribution_no_approx = clmm.theory.compute_reduced_tangential_shear(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            model_z_distrib_dict['func'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='distribution', beta_kwargs={'zmax':zsrc_max}, approx=None)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### **Cases 3 : Redshift distribution and approximation**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### If we want a faster approach we can use an approximation for the reduced shear, using 1 or 2 order of Taylor expansion for the expression of the reduced shear."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "gt_distribution_order1 = clmm.theory.compute_reduced_tangential_shear(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            model_z_distrib_dict['func'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='distribution', beta_kwargs={'zmax':zsrc_max},  approx=\"order1\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "gt_distribution_order2 = clmm.theory.compute_reduced_tangential_shear(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            model_z_distrib_dict['func'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='distribution', beta_kwargs={'zmax':zsrc_max},  approx=\"order2\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### **Cases 4 : Mean lensing efficiencies and approximation**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally, we can also model the reduced shear if the only information we have about the source redshift distribution is through the mean lensing efficiency parameters $\\langle\\beta_s\\rangle$ and $\\langle\\beta_s^2\\rangle$. \\\n",
+    "In this case, we need to use an approximation for the formula. This is the fastest approach."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "gt_beta_1 = clmm.theory.compute_reduced_tangential_shear(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', beta_kwargs={'zmax':zsrc_max}, z_src_info='beta', approx=\"order1\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "gt_beta_2 = clmm.theory.compute_reduced_tangential_shear(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', beta_kwargs={'zmax':zsrc_max}, z_src_info='beta', approx=\"order2\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 3 - Comparison of the four cases for the reduced tangnetial shear"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, (ax1, ax2) = plt.subplots(1,2,sharex=True)\n",
+    "\n",
+    "ax1.loglog(gc_object.profile['radius'], gc_object.profile['gt'], 'gs-', label = 'data')\n",
+    "\n",
+    "ax1.loglog(rr, gt_discrete, 'k.-', label = 'discrete')\n",
+    "ax1.loglog(rr, gt_distribution_no_approx, 'r.-', label = 'distribution, no approx')\n",
+    "ax1.loglog(rr, gt_distribution_order1, 'ro-', label = 'distribution, order 1 approx')\n",
+    "ax1.loglog(rr, gt_distribution_order2, 'rd-', label = 'distribution, order 2 approx')\n",
+    "ax1.loglog(rr, gt_beta_1, 'b--', label = 'beta, order 1 approx')\n",
+    "ax1.loglog(rr, gt_beta_2, 'bx-', label = 'beta, order 2 approx')\n",
+    "\n",
+    "ax1.set_ylabel('$g_t$')\n",
+    "ax1.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax1.legend()\n",
+    "\n",
+    "ax2.plot(rr, 100*(gt_discrete-gt_discrete)/gt_discrete, 'k.-', label = 'discrete')\n",
+    "ax2.plot(rr, 100*(gt_distribution_no_approx-gt_discrete)/gt_discrete, 'r-.', label = 'distribution, no approx')\n",
+    "ax2.plot(rr, 100*(gt_distribution_order1-gt_discrete)/gt_discrete, 'ro-', label = 'distribution, order 1 approx')\n",
+    "ax2.plot(rr, 100*(gt_distribution_order2-gt_discrete)/gt_discrete, 'rd-', label = 'distribution, order 2 approx')\n",
+    "ax2.plot(rr, 100*(gt_beta_1-gt_discrete)/gt_discrete, 'b--', label = 'beta, order 1 approx')\n",
+    "ax2.plot(rr, 100*(gt_beta_2-gt_discrete)/gt_discrete, 'bx-', label = 'beta, order 2 approx')\n",
+    "\n",
+    "ax2.set_ylabel('%')\n",
+    "ax2.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax2.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "All modeled profiles give similar results. They do not correspond to the profile computed from the data in the inner part, because of different ways of constructing the profiles (taking the average radial point and reduced shear value in a bin or computing the average expected reduced shear at a given radius).\n",
+    "\n",
+    "The profiles computed using an approximation for the reduced shear formula are lower by a few percents, especially in the inner region. The profiles computed from a redshift distribution or a known source redshift differ at the subpercent level."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 4- Modeling and plotting the other quantities"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here we will just compute models for cases 1, 2 and 4. For the shear and convergence there is no need for an approximated formula. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gammat_discrete = np.zeros(rr.size)\n",
+    "\n",
+    "for i in range(rr.size):\n",
+    "    gammat_discrete[i] = np.mean(clmm.theory.compute_tangential_shear(rr[i], cluster_m, concentration, cluster_z, \n",
+    "                                            source_catalog['z'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='discrete'))\n",
+    "    \n",
+    "gammat_distribution = clmm.theory.compute_tangential_shear(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            model_z_distrib_dict['func'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='distribution', beta_kwargs={'zmax':zsrc_max})\n",
+    "    \n",
+    "gammat_beta = clmm.theory.compute_tangential_shear(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', z_src_info='beta', beta_kwargs={'zmax':zsrc_max})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, (ax1, ax2) = plt.subplots(1,2,sharex=True)\n",
+    "\n",
+    "ax1.loglog(rr, gammat_discrete, 'k.-', label = 'discrete')\n",
+    "ax1.loglog(rr, gammat_distribution, 'rx-', label = 'distribution, no approx')\n",
+    "ax1.loglog(rr, gammat_beta, 'b--', label = 'beta, no approx')\n",
+    "\n",
+    "ax1.set_ylabel('$\\\\gamma_t$')\n",
+    "ax1.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax1.legend()\n",
+    "\n",
+    "ax2.plot(rr, 100*(gammat_discrete-gammat_discrete)/gammat_discrete, 'k.-', label = 'discrete')\n",
+    "ax2.plot(rr, 100*(gammat_distribution-gammat_discrete)/gammat_discrete, 'rx-', label = 'distribution, no approx')\n",
+    "ax2.plot(rr, 100*(gammat_beta-gammat_discrete)/gammat_discrete, 'b--', label = 'beta, no approx')\n",
+    "\n",
+    "ax2.set_ylabel('%')\n",
+    "ax2.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax2.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "kappa_discrete = np.zeros(rr.size)\n",
+    "\n",
+    "for i in range(rr.size):\n",
+    "    kappa_discrete[i] = np.mean(clmm.theory.compute_convergence(rr[i], cluster_m, concentration, cluster_z, \n",
+    "                                            source_catalog['z'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='discrete'))\n",
+    "    \n",
+    "kappa_distribution = clmm.theory.compute_convergence(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            model_z_distrib_dict['func'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='distribution', beta_kwargs={'zmax':zsrc_max})\n",
+    "    \n",
+    "kappa_beta = clmm.theory.compute_convergence(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', z_src_info='beta', beta_kwargs={'zmax':zsrc_max})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, (ax1, ax2) = plt.subplots(1,2,sharex=True)\n",
+    "\n",
+    "ax1.loglog(rr, kappa_discrete, 'k.-', label = 'discrete')\n",
+    "ax1.loglog(rr, kappa_distribution, 'rx-', label = 'distribution, no approx')\n",
+    "ax1.loglog(rr, kappa_beta, 'b--', label = 'beta, no approx')\n",
+    "\n",
+    "ax1.set_ylabel('$\\\\kappa$')\n",
+    "ax1.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax1.legend()\n",
+    "\n",
+    "ax2.plot(rr, 100*(kappa_discrete-kappa_discrete)/kappa_discrete, 'k.-', label = 'discrete')\n",
+    "ax2.plot(rr, 100*(kappa_distribution-kappa_discrete)/kappa_discrete, 'rx-', label = 'distribution, no approx')\n",
+    "ax2.plot(rr, 100*(kappa_beta-kappa_discrete)/kappa_discrete, 'b--', label = 'beta, no approx')\n",
+    "\n",
+    "ax2.set_ylabel('%')\n",
+    "ax2.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax2.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mu_discrete = np.zeros(rr.size)\n",
+    "\n",
+    "for i in range(rr.size):\n",
+    "    mu_discrete[i] = np.mean(clmm.theory.compute_magnification(rr[i], cluster_m, concentration, cluster_z, \n",
+    "                                            source_catalog['z'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='discrete'))\n",
+    "    \n",
+    "mu_distribution = clmm.theory.compute_magnification(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            model_z_distrib_dict['func'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='distribution', approx=None, beta_kwargs={'zmax':zsrc_max})\n",
+    "    \n",
+    "mu_beta_1 = clmm.theory.compute_magnification(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', z_src_info='beta',  approx=\"order1\", beta_kwargs={'zmax':zsrc_max})\n",
+    "\n",
+    "mu_beta_2 = clmm.theory.compute_magnification(rr, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', z_src_info='beta',  approx=\"order2\", beta_kwargs={'zmax':zsrc_max})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, (ax1, ax2) = plt.subplots(1,2,sharex=True)\n",
+    "\n",
+    "ax1.loglog(rr, mu_discrete, 'k.-', label = 'discrete')\n",
+    "ax1.loglog(rr, mu_distribution, 'rs--', label = 'distribution, no approx')\n",
+    "ax1.loglog(rr, mu_beta_1, 'bx-', label = 'beta, order 1 approx')\n",
+    "ax1.loglog(rr, mu_beta_2, 'b--', label = 'beta, order 2 approx')\n",
+    "\n",
+    "ax1.set_ylabel('$\\\\mu$')\n",
+    "ax1.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax1.legend()\n",
+    "\n",
+    "ax2.plot(rr, 100*(mu_discrete-mu_discrete)/mu_discrete, 'k.-', label = 'discrete')\n",
+    "ax2.plot(rr, 100*(mu_distribution-mu_discrete)/mu_discrete, 'rs--', label = 'distribution, no approx')\n",
+    "ax2.plot(rr, 100*(mu_beta_1-mu_discrete)/mu_discrete, 'bx-', label = 'beta, order 1 approx')\n",
+    "ax2.plot(rr, 100*(mu_beta_2-mu_discrete)/mu_discrete, 'b--', label = 'beta, order 2 approx')\n",
+    "\n",
+    "ax2.set_ylabel('%')\n",
+    "ax2.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax2.legend()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mu_bias_discrete = np.zeros(rr.size)\n",
+    "\n",
+    "alpha = 2.7 #arbitrary value for the slope of the number density\n",
+    "\n",
+    "for i in range(rr.size):\n",
+    "    mu_bias_discrete[i] = np.mean(clmm.theory.compute_magnification_bias(rr[i], alpha, cluster_m, concentration, cluster_z, \n",
+    "                                            source_catalog['z'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='discrete'))\n",
+    "    \n",
+    "mu_bias_distribution = clmm.theory.compute_magnification_bias(rr, alpha, cluster_m, concentration, cluster_z, \n",
+    "                                            model_z_distrib_dict['func'], cosmo, delta_mdef=500,  massdef='critical', z_src_info='distribution', approx=None, beta_kwargs={'zmax':zsrc_max})\n",
+    "    \n",
+    "mu_bias_beta_1 = clmm.theory.compute_magnification_bias(rr, alpha, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', z_src_info='beta',  approx=\"order1\", beta_kwargs={'zmax':zsrc_max})\n",
+    "\n",
+    "mu_bias_beta_2 = clmm.theory.compute_magnification_bias(rr, alpha, cluster_m, concentration, cluster_z, \n",
+    "                                            [beta_s_mean, beta_s_square_mean], cosmo, delta_mdef=500,  massdef='critical', z_src_info='beta',  approx=\"order2\", beta_kwargs={'zmax':zsrc_max})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, (ax1, ax2) = plt.subplots(1,2,sharex=True)\n",
+    "\n",
+    "ax1.loglog(rr, mu_bias_discrete, 'k.-', label = 'discrete')\n",
+    "ax1.loglog(rr, mu_bias_distribution, 'rs--', label = 'distribution, no approx')\n",
+    "ax1.loglog(rr, mu_bias_beta_1, 'bx-', label = 'beta, order 1 approx')\n",
+    "ax1.loglog(rr, mu_bias_beta_2, 'b--', label = 'beta, order 2 approx')\n",
+    "\n",
+    "\n",
+    "ax1.set_ylabel('$\\\\delta_{\\\\mu} + 1$')\n",
+    "ax1.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax1.legend()\n",
+    "\n",
+    "ax2.plot(rr, 100*(mu_bias_discrete-mu_bias_discrete)/mu_bias_discrete, 'k.-', label = 'discrete')\n",
+    "ax2.plot(rr, 100*(mu_bias_distribution-mu_bias_discrete)/mu_bias_discrete, 'rs--', label = 'distribution, no approx')\n",
+    "ax2.plot(rr, 100*(mu_bias_beta_1-mu_bias_discrete)/mu_bias_discrete, 'bx-', label = 'beta, order 1 approx')\n",
+    "ax2.plot(rr, 100*(mu_bias_beta_2-mu_bias_discrete)/mu_bias_discrete, 'b--', label = 'beta, order 2 approx')\n",
+    "\n",
+    "ax2.set_ylabel('%')\n",
+    "ax2.set_xlabel('radius [Mpc]')\n",
+    "\n",
+    "ax2.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**When making approximation, there may be large differences (~few 10% level) in the inner regions (compare to when using the exact redshift of the sources), especially for magnification and magnification bias. However, these approaches are very fast to compute. The user as to chose the appropriate method depending on the use case.**"
+   ]
+  }
+ ],
+ "metadata": {
+  "anaconda-cloud": {},
+  "kernelspec": {
+   "display_name": "wrk",
+   "language": "python",
+   "name": "wrk"
+  },
+  "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.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/demo_theory_functionality.ipynb b/examples/demo_theory_functionality.ipynb
index ae0a707ef..35259f27e 100644
--- a/examples/demo_theory_functionality.ipynb
+++ b/examples/demo_theory_functionality.ipynb
@@ -103,7 +103,7 @@
     "\n",
     "#source properties\n",
     "z_source = 2. #all sources in the same plan\n",
-    "z_distrib_func = clmm.utils._chang_z_distrib #sources redshift following a distribution\n",
+    "z_distrib_func = clmm.z_distributions.chang2013 #sources redshift following a distribution\n",
     "alpha = [2, -0.5] "
    ]
   },
@@ -173,7 +173,7 @@
     "gammat = m.compute_tangential_shear(r3d, mdelta=cluster_mass, cdelta=cluster_concentration, z_cluster=z_cl, \n",
     "                                    z_source=z_source, cosmo=cosmo, delta_mdef=mass_Delta, \n",
     "                                    halo_profile_model=density_profile_parametrization, \n",
-    "                                    z_src_model='single_plane')"
+    "                                    z_src_info='discrete')"
    ]
   },
   {
@@ -186,7 +186,7 @@
     "                               z_cluster=z_cl, z_source=z_source,\n",
     "                               cosmo=cosmo, delta_mdef=mass_Delta, \n",
     "                               halo_profile_model=density_profile_parametrization, \n",
-    "                               z_src_model='single_plane')"
+    "                               z_src_info='discrete')"
    ]
   },
   {
@@ -199,8 +199,7 @@
     "                                        z_cluster=z_cl, z_source=z_source, cosmo=cosmo, \n",
     "                                        delta_mdef=mass_Delta, \n",
     "                                        halo_profile_model=density_profile_parametrization, \n",
-    "                                        z_src_model='single_plane')\n",
-    "        "
+    "                                        z_src_info='discrete')"
    ]
   },
   {
@@ -213,7 +212,7 @@
     "                                        z_cluster=z_cl, z_source=z_source, cosmo=cosmo, \n",
     "                                        delta_mdef=mass_Delta, \n",
     "                                        halo_profile_model=density_profile_parametrization, \n",
-    "                                        z_src_model='single_plane')\n"
+    "                                        z_src_info='discrete')"
    ]
   },
   {
@@ -226,7 +225,7 @@
     "                                        z_cluster=z_cl, z_source=z_source, cosmo=cosmo, \n",
     "                                        delta_mdef=mass_Delta, \n",
     "                                        halo_profile_model=density_profile_parametrization, \n",
-    "                                        z_src_model='single_plane')\n"
+    "                                        z_src_info='discrete')\n"
    ]
   },
   {
@@ -243,11 +242,11 @@
    "outputs": [],
    "source": [
     "gt_z = m.compute_reduced_tangential_shear(r3d, mdelta=cluster_mass, cdelta=cluster_concentration, \n",
-    "                                        z_cluster=z_cl, z_source=z_source, cosmo=cosmo, \n",
+    "                                        z_cluster=z_cl, z_source=z_distrib_func, cosmo=cosmo, \n",
     "                                        delta_mdef=mass_Delta, \n",
     "                                        halo_profile_model=density_profile_parametrization,\n",
-    "                                        z_src_model='schrabback18',\n",
-    "                                        z_distrib_func = z_distrib_func)"
+    "                                        z_src_info='distribution',\n",
+    "                                        approx='order2')"
    ]
   },
   {
@@ -422,7 +421,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.10.9"
   }
  },
  "nbformat": 4,
diff --git a/examples/demo_theory_functionality_oo.ipynb b/examples/demo_theory_functionality_oo.ipynb
index d1c816267..07f79a7df 100644
--- a/examples/demo_theory_functionality_oo.ipynb
+++ b/examples/demo_theory_functionality_oo.ipynb
@@ -103,7 +103,7 @@
     "\n",
     "#source properties\n",
     "z_source = 2. #all sources in the same plane\n",
-    "z_distrib_func = clmm.utils._chang_z_distrib #sources redshift following a distribution\n",
+    "z_distrib_func = clmm.z_distributions.chang2013 #sources redshift following a distribution\n",
     "alpha = [2, -0.5] "
    ]
   },
@@ -130,8 +130,8 @@
     "\n",
     "gt = moo.eval_reduced_tangential_shear(r3d, z_cl, z_source)\n",
     "#Lensing quantities assuming sources follow a given redshift distribution.\n",
-    "gt_z = moo.eval_reduced_tangential_shear(r3d, z_cl, z_source, z_src_model='schrabback18',\n",
-    "                                        z_distrib_func = z_distrib_func)\n",
+    "gt_z = moo.eval_reduced_tangential_shear(r3d, z_cl, z_distrib_func, z_src_info='distribution',\n",
+    "                                        approx='order2')\n",
     "\n",
     "mu = moo.eval_magnification(r3d, z_cl, z_source)\n",
     "mu_bias = moo.eval_magnification_bias(r3d, z_cl, z_source, alpha)"
@@ -293,9 +293,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "wrk",
    "language": "python",
-   "name": "python3"
+   "name": "wrk"
   },
   "language_info": {
    "codemirror_mode": {
@@ -307,7 +307,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.10"
+   "version": "3.10.6"
   }
  },
  "nbformat": 4,


Test patch:
diff --git a/tests/test_mockdata.py b/tests/test_mockdata.py
index 8e988e076..64ad8e50e 100644
--- a/tests/test_mockdata.py
+++ b/tests/test_mockdata.py
@@ -6,7 +6,7 @@
 import clmm.dataops as da
 from clmm.support import mock_data as mock
 from clmm.support.sampler import fitters
-from clmm.utils import _chang_z_distrib, _srd_z_distrib
+from clmm import z_distributions as zdist
 
 TOLERANCE = {'rtol': 5.0e-4, 'atol': 1.e-4}
 
@@ -148,9 +148,9 @@ def test_z_distr():
     assert_equal(np.count_nonzero((data['z'] < zmin) | (data['z'] > zmax)), 0)
     # Check that the z distribution follows Chang13 distribution
     hist = np.histogram(data['z'], bins=bins)
-    norm = _chang_z_distrib(zmax, is_cdf=True)-_chang_z_distrib(zmin, is_cdf=True)
+    norm = zdist.chang2013(zmax, is_cdf=True)-zdist.chang2013(zmin, is_cdf=True)
     # Expected number of galaxies in bin i = Ntot*distrib(z_{bin center})*binsize/norm
-    chang = np.array([mock._chang_z_distrib(z)*ngals *
+    chang = np.array([zdist.chang2013(z)*ngals *
                      0.1/norm for z in bins[:-1]+0.05])
     assert_allclose(hist[0], chang, atol=100, rtol=0.1)
 
@@ -161,9 +161,9 @@ def test_z_distr():
     assert_equal(np.count_nonzero((data['z'] < zmin) | (data['z'] > zmax)), 0)
     # Check that the z distribution follows Chang13 distribution
     hist = np.histogram(data['z'], bins=bins)
-    norm = _srd_z_distrib(zmax, is_cdf=True)-_srd_z_distrib(zmin, is_cdf=True)
+    norm = zdist.desc_srd(zmax, is_cdf=True)-zdist.desc_srd(zmin, is_cdf=True)
     # Expected number of galaxies in bin i = Ntot*distrib(z_{bin center})*binsize/norm
-    srd = np.array([mock._srd_z_distrib(z)*ngals*0.1/norm for z in bins[:-1]+0.05])
+    srd = np.array([zdist.desc_srd(z)*ngals*0.1/norm for z in bins[:-1]+0.05])
     assert_allclose(hist[0],srd,atol=100,rtol=0.1)
 
 
diff --git a/tests/test_theory.py b/tests/test_theory.py
index 78fb0c76f..067cc83d6 100644
--- a/tests/test_theory.py
+++ b/tests/test_theory.py
@@ -7,6 +7,8 @@
 from clmm.constants import Constants as clc
 from clmm.galaxycluster import GalaxyCluster
 from clmm import GCData
+from clmm.utils import compute_beta_s_square_mean, compute_beta_s_mean
+from clmm.z_distributions import chang2013, desc_srd
 
 TOLERANCE = {'rtol': 1.0e-8}
 
@@ -94,7 +96,7 @@ def load_validation_config(halo_profile_model=None):
         'z_source': testcase['z_source'],
         'delta_mdef': testcase['mass_Delta'],
         'halo_profile_model': testcase['density_profile_parametrization'],
-        'z_src_model': 'single_plane',
+        'z_src_info': 'discrete',
     }
 
     return {'TEST_CASE': testcase, 'z_source': testcase['z_source'],
@@ -288,16 +290,13 @@ def test_profiles(modeling_data, profile_init):
             alpha_ein = None
 
         assert_allclose(
-            mod.eval_3d_density(cfg['RHO_PARAMS']['r3d'],
-                                cfg['RHO_PARAMS']['z_cl'], verbose=True),
+            mod.eval_3d_density(cfg['RHO_PARAMS']['r3d'], cfg['RHO_PARAMS']['z_cl'], verbose=True),
             cfg['numcosmo_profiles']['rho'], reltol)
         assert_allclose(
-            mod.eval_surface_density(cfg['SIGMA_PARAMS']['r_proj'],
-                                     cfg['SIGMA_PARAMS']['z_cl'], verbose=True),
+            mod.eval_surface_density(cfg['SIGMA_PARAMS']['r_proj'], cfg['SIGMA_PARAMS']['z_cl'], verbose=True),
             cfg['numcosmo_profiles']['Sigma'], reltol)
         assert_allclose(
-            mod.eval_excess_surface_density(cfg['SIGMA_PARAMS']['r_proj'],
-                                            cfg['SIGMA_PARAMS']['z_cl'], verbose=True),
+            mod.eval_excess_surface_density(cfg['SIGMA_PARAMS']['r_proj'], cfg['SIGMA_PARAMS']['z_cl'], verbose=True),
             cfg['numcosmo_profiles']['DeltaSigma'], reltol)
         if mod.backend == 'ct':
             assert_raises(ValueError, mod.eval_excess_surface_density,
@@ -305,14 +304,11 @@ def test_profiles(modeling_data, profile_init):
 
         # Functional interface tests
         # alpha_ein is None unless testing Einasto with the NC backend
-        assert_allclose(theo.compute_3d_density(cosmo=cosmo, **cfg['RHO_PARAMS'],
-                                                alpha_ein=alpha_ein, verbose=True),
+        assert_allclose(theo.compute_3d_density(cosmo=cosmo, **cfg['RHO_PARAMS'], alpha_ein=alpha_ein, verbose=True),
                         cfg['numcosmo_profiles']['rho'], reltol)
-        assert_allclose(theo.compute_surface_density(cosmo=cosmo, **cfg['SIGMA_PARAMS'],
-                                                     alpha_ein=alpha_ein, verbose=True),
+        assert_allclose(theo.compute_surface_density(cosmo=cosmo, **cfg['SIGMA_PARAMS'], alpha_ein=alpha_ein, verbose=True),
                         cfg['numcosmo_profiles']['Sigma'], reltol)
-        assert_allclose(theo.compute_excess_surface_density(cosmo=cosmo, **cfg['SIGMA_PARAMS'],
-                                                            alpha_ein=alpha_ein, verbose=True),
+        assert_allclose(theo.compute_excess_surface_density(cosmo=cosmo, **cfg['SIGMA_PARAMS'], alpha_ein=alpha_ein, verbose=True),
                         cfg['numcosmo_profiles']['DeltaSigma'], reltol)
 
     else:
@@ -433,17 +429,16 @@ def helper_physics_functions(func, additional_kwargs={}):
     kwargs.update(additional_kwargs)
 
     # Test defaults
-
     defaulttruth = func(**kwargs, delta_mdef=200, halo_profile_model='nfw',
-                        z_src_model='single_plane')
+                        z_src_info='discrete')
     assert_allclose(func(**kwargs, delta_mdef=200), defaulttruth, **TOLERANCE)
     assert_allclose(func(**kwargs, halo_profile_model='nfw'), defaulttruth, **TOLERANCE)
-    assert_allclose(func(**kwargs, z_src_model='single_plane'), defaulttruth, **TOLERANCE)
+    assert_allclose(func(**kwargs, z_src_info='discrete'), defaulttruth, **TOLERANCE)
 
-    # Test for exception on unsupported z_src_model and halo profiles
+    # Test for exception on unsupported z_src_info and halo profiles
     assert_raises(ValueError, func, **kwargs, halo_profile_model='blah')
     assert_raises(ValueError, func, **kwargs, massdef='blah')
-    assert_raises(ValueError, func, **kwargs, z_src_model='blah')
+    assert_raises(ValueError, func, **kwargs, z_src_info='blah')
 
 
 def test_shear_convergence_unittests(modeling_data, profile_init):
@@ -492,64 +487,228 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                               theo.compute_magnification, cosmo=cosmo, **cfg['GAMMA_PARAMS'])
                 del cfg['GAMMA_PARAMS']['alpha_ein']
 
-        # Validate tangential shear
+        # Validate tangential shear - discrete case
         gammat = theo.compute_tangential_shear(cosmo=cosmo, **cfg['GAMMA_PARAMS'])
         assert_allclose(gammat, cfg['numcosmo_profiles']['gammat'], reltol)
 
-        # Validate convergence
+        # Validate convergence - discrete case
         kappa = theo.compute_convergence(cosmo=cosmo, **cfg['GAMMA_PARAMS'])
         assert_allclose(kappa, cfg['numcosmo_profiles']['kappa'], reltol)
 
-        # Validate reduced tangential shear
+        # Validate reduced tangential shear - discrete case
+        gt = theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg['GAMMA_PARAMS'])
+        assert_allclose(gt, gammat/(1.0-kappa), 1.0e-10)
+        assert_allclose(gt, cfg['numcosmo_profiles']['gt'], 1.e2*reltol)
+
+        # Validate magnification - discrete case
+        mu = theo.compute_magnification(cosmo=cosmo, **cfg['GAMMA_PARAMS'])
+        assert_allclose(mu, 1./((1-kappa)**2-abs(gammat)**2), 1.0e-10)
+        assert_allclose(mu, cfg['numcosmo_profiles']['mu'], 1.e2*reltol)
+
+        # Validate magnification bias - discrete case
+        # this if will be removed once compute_magnification_bias takes alpha_ein
+        if profile_init=='einasto' and theo.be_nick=='nc':
+            del cfg['GAMMA_PARAMS']['alpha_ein']
+            cfg['GAMMA_PARAMS']['verbose'] = True
+
+        alpha = 3.78
         assert_allclose(
-            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg['GAMMA_PARAMS']),
-            gammat/(1.0-kappa), 1.0e-10)
+            theo.compute_magnification_bias(cosmo=cosmo, **cfg['GAMMA_PARAMS'], alpha=alpha),
+            (1./((1-kappa)**2-abs(gammat)**2))**(alpha - 1), 1.0e-10)
         assert_allclose(
-            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg['GAMMA_PARAMS']),
-            cfg['numcosmo_profiles']['gt'], 1.e2*reltol)
+            theo.compute_magnification_bias(cosmo=cosmo, **cfg['GAMMA_PARAMS'], alpha=alpha),
+            cfg['numcosmo_profiles']['mu']**(alpha - 1), 1.e3*reltol)
 
-        beta_s_mean, beta_s_square_mean = 0.9, 0.6
         cfg_inf = load_validation_config()
+
+        # compute some values
         cfg_inf['GAMMA_PARAMS']['z_source'] = 1000.
+        beta_s_mean = compute_beta_s_mean(
+            cfg_inf['GAMMA_PARAMS']['z_cluster'], cfg_inf['GAMMA_PARAMS']['z_source'], cosmo)
+        beta_s_square_mean = compute_beta_s_square_mean(
+            cfg_inf['GAMMA_PARAMS']['z_cluster'], cfg_inf['GAMMA_PARAMS']['z_source'], cosmo)
+
         gammat_inf = theo.compute_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])
         kappa_inf = theo.compute_convergence(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])
 
-        cfg_inf['GAMMA_PARAMS']['z_src_model'] = 'applegate14'
+        # test z_src = chang2013 distribution
+        cfg_inf['GAMMA_PARAMS']['z_source'] = chang2013
+        cfg_inf['GAMMA_PARAMS']['z_src_info'] = 'distribution'
+
+        # store original values
+        r_proj = cfg_inf['GAMMA_PARAMS']['r_proj']
+        # use only 5 largest radii
+        cfg_inf['GAMMA_PARAMS']['r_proj'] = cfg_inf['GAMMA_PARAMS']['r_proj'][-5:]
+        # calculate some true values
+        gt = theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])
+        mu = theo.compute_magnification(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])
+        mu_bias = theo.compute_magnification_bias(
+            cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'], alpha=alpha)
+        cfg_inf['GAMMA_PARAMS']['r_proj'] = r_proj
+        # tangential shear
+        assert_allclose(
+            theo.compute_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            beta_s_mean * gammat_inf, 1.0e-10)
+
+        # convergence
+        assert_allclose(
+            theo.compute_convergence(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            beta_s_mean * kappa_inf, 1.0e-10)
+
+        # reduced tangential shear
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order1'
         assert_allclose(
-            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
-                                                  beta_s_mean=beta_s_mean,
-                                                  beta_s_square_mean=beta_s_square_mean),
-            beta_s_mean * gammat_inf/(1.0 - beta_s_square_mean / beta_s_mean * kappa_inf),
+            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            beta_s_mean * gammat_inf/(1.0 - beta_s_mean * kappa_inf),
             1.0e-10)
+        assert_allclose(
+            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])[-5:],
+            gt, 3.0e-6)
 
-        cfg_inf['GAMMA_PARAMS']['z_src_model'] = 'schrabback18'
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order2'
         assert_allclose(
-            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
-                                                  beta_s_mean=beta_s_mean,
-                                                  beta_s_square_mean=beta_s_square_mean),
+            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
             (1. + (beta_s_square_mean / (beta_s_mean * beta_s_mean) - 1.) * beta_s_mean *
              kappa_inf) * (beta_s_mean * gammat_inf / (1. - beta_s_mean * kappa_inf)),
             1.0e-10)
+        assert_allclose(
+            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])[-5:],
+            gt, 2.0e-6)
 
-        # Validate magnification
-        assert_allclose(theo.compute_magnification(cosmo=cosmo, **cfg['GAMMA_PARAMS']),
-                        1./((1-kappa)**2-abs(gammat)**2), 1.0e-10)
-        assert_allclose(theo.compute_magnification(cosmo=cosmo, **cfg['GAMMA_PARAMS']),
-                        cfg['numcosmo_profiles']['mu'], 1.e2*reltol)
+        # magnification
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order1'
+        assert_allclose(
+            theo.compute_magnification(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            1+2*beta_s_mean*kappa_inf, 1.0e-10)
+        assert_allclose(
+            theo.compute_magnification(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])[-5:],
+            mu, 2.0e-8)
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order2'
+        assert_allclose(
+            theo.compute_magnification(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            1+2*beta_s_mean*kappa_inf+beta_s_square_mean*gammat_inf**2\
+            +3*beta_s_square_mean*kappa_inf**2, 1.0e-10)
+        assert_allclose(
+            theo.compute_magnification(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])[-5:],
+            mu, 1.0e-10)
 
-        # Validate magnification bias
-        alpha = 3.78
+        # magnification bias
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order1'
+        assert_allclose(
+            theo.compute_magnification_bias(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'], alpha=alpha),
+            1+(alpha-1)*(2*beta_s_mean*kappa_inf), 1.0e-10)
+        assert_allclose(
+            theo.compute_magnification_bias(
+                cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'], alpha=alpha)[-5:],
+            mu_bias, 4.0e-8)
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order2'
+        assert_allclose(
+            theo.compute_magnification_bias(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'], alpha=alpha),
+            1+(alpha-1)*(2*beta_s_mean*kappa_inf+beta_s_square_mean*gammat_inf**2)\
+            +(2*alpha-1)*(alpha-1)*beta_s_square_mean*kappa_inf**2, 1.0e-10)
+        assert_allclose(
+            theo.compute_magnification_bias(
+                cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'], alpha=alpha)[-5:],
+            mu_bias, 1.0e-10)
+
+        # test errors and also prepare for the next round of tests
+        del cfg_inf['GAMMA_PARAMS']['approx']
+        # test ValueError from unsupported approx
+        assert_raises(ValueError,
+                      theo.compute_reduced_tangential_shear, cosmo=cosmo,
+                      **cfg_inf['GAMMA_PARAMS'],
+                      approx='notvalid')
+        assert_raises(ValueError,
+                      theo.compute_magnification, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      approx='notvalid')
+        assert_raises(ValueError,
+                      theo.compute_magnification_bias, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      alpha=alpha,
+                      approx='notvalid')
+        # test KeyError from invalid key in beta_kwargs
+        assert_raises(KeyError,
+                      theo.compute_tangential_shear, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      beta_kwargs={'notavalidkey': 0.})
+        assert_raises(KeyError,
+                      theo.compute_convergence, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      beta_kwargs={'notavalidkey': 0.})
+        assert_raises(KeyError,
+                      theo.compute_reduced_tangential_shear, cosmo=cosmo,
+                      **cfg_inf['GAMMA_PARAMS'],
+                      beta_kwargs={'notavalidkey': 0.})
+        assert_raises(KeyError,
+                      theo.compute_magnification, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      beta_kwargs={'notavalidkey': 0.})
+        assert_raises(KeyError,
+                      theo.compute_magnification_bias, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      alpha=alpha,
+                      beta_kwargs={'notavalidkey': 0.})
+        # test ValueError from unsupported z_src_info
+        cfg_inf['GAMMA_PARAMS']['z_src_info'] = 'notvalid'
+        assert_raises(ValueError,
+                      theo.compute_tangential_shear, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])
+        assert_raises(ValueError,
+                      theo.compute_convergence, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'])
+        assert_raises(ValueError,
+                      theo.compute_reduced_tangential_shear, cosmo=cosmo,
+                      **cfg_inf['GAMMA_PARAMS'],
+                      approx='order1')
+        assert_raises(ValueError,
+                      theo.compute_magnification, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      approx='order1')
+        assert_raises(ValueError,
+                      theo.compute_magnification_bias, cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'],
+                      alpha=2,
+                      approx='order1')
+
+        #test z_src_info = 'beta'
+        beta_s_mean, beta_s_square_mean = 0.9, 0.6
+        cfg_inf['GAMMA_PARAMS']['z_source'] = (beta_s_mean, beta_s_square_mean)
+        cfg_inf['GAMMA_PARAMS']['z_src_info'] = 'beta'
+        # tangential shear
+        assert_allclose(
+            theo.compute_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            beta_s_mean * gammat_inf, 1.0e-10)
 
-        # this if will be removed once compute_magnification_bias takes alpha_ein
-        if profile_init=='einasto' and theo.be_nick=='nc':
-            del cfg['GAMMA_PARAMS']['alpha_ein']
+        # convergence
+        assert_allclose(
+            theo.compute_convergence(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            beta_s_mean * kappa_inf, 1.0e-10)
 
+        # reduced tangential shear
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order1'
         assert_allclose(
-            theo.compute_magnification_bias(cosmo=cosmo, **cfg['GAMMA_PARAMS'], alpha=alpha),
-            (1./((1-kappa)**2-abs(gammat)**2))**(alpha - 1), 1.0e-10)
+            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            beta_s_mean * gammat_inf/(1.0 - beta_s_mean * kappa_inf),
+            1.0e-10)
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order2'
         assert_allclose(
-            theo.compute_magnification_bias(cosmo=cosmo, **cfg['GAMMA_PARAMS'], alpha=alpha),
-            cfg['numcosmo_profiles']['mu']**(alpha - 1), 1.e3*reltol)
+            theo.compute_reduced_tangential_shear(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            (1. + (beta_s_square_mean / (beta_s_mean * beta_s_mean) - 1.) * beta_s_mean *
+             kappa_inf) * (beta_s_mean * gammat_inf / (1. - beta_s_mean * kappa_inf)),
+            1.0e-10)
+
+        # magnification
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order1'
+        assert_allclose(
+            theo.compute_magnification(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            1+2*beta_s_mean*kappa_inf, 1.0e-10)
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order2'
+        assert_allclose(
+            theo.compute_magnification(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS']),
+            1+2*beta_s_mean*kappa_inf+beta_s_square_mean*gammat_inf**2\
+            +3*beta_s_square_mean*kappa_inf**2, 1.0e-10)
+
+        # magnification bias
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order1'
+        assert_allclose(
+            theo.compute_magnification_bias(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'], alpha=alpha),
+            1+(alpha-1)*(2*beta_s_mean*kappa_inf), 1.0e-10)
+        cfg_inf['GAMMA_PARAMS']['approx'] = 'order2'
+        assert_allclose(
+            theo.compute_magnification_bias(cosmo=cosmo, **cfg_inf['GAMMA_PARAMS'], alpha=alpha),
+            1+(alpha-1)*(2*beta_s_mean*kappa_inf+beta_s_square_mean*gammat_inf**2)\
+            +(2*alpha-1)*(alpha-1)*beta_s_square_mean*kappa_inf**2, 1.0e-10)
 
         # Check that shear, reduced shear and convergence return zero
         # and magnification and magnification bias return one
@@ -627,9 +786,9 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
         if profile_init=='einasto' and theo.be_nick=='nc':
             mod.set_einasto_alpha(cfg['TEST_CASE']['alpha_einasto'])
 
-        # Validate tangential shear
         profile_pars = [cfg['GAMMA_PARAMS']['r_proj'], cfg['GAMMA_PARAMS']['z_cluster'],
                         cfg['GAMMA_PARAMS']['z_source']]
+        # Validate tangential shear
         gammat = mod.eval_tangential_shear(*profile_pars)
         assert_allclose(gammat,
                         cfg['numcosmo_profiles']['gammat'], reltol)
@@ -645,36 +804,76 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
         assert_allclose(mod.eval_reduced_tangential_shear(*profile_pars),
                         cfg['numcosmo_profiles']['gt'], 1.e2*reltol)
 
+        # Validate magnification
+        assert_allclose(mod.eval_magnification(*profile_pars),
+                        1./((1-kappa)**2-abs(gammat)**2), 1.0e-10)
+        assert_allclose(mod.eval_magnification(*profile_pars),
+                        cfg['numcosmo_profiles']['mu'], 1.e2*reltol)
+
+        # Validate magnification bias
+        alpha = -1.78
+        assert_allclose(mod.eval_magnification_bias(*profile_pars, alpha=alpha),
+                        1./((1-kappa)**2-abs(gammat)**2)**(alpha-1), 1.0e-10)
+        assert_allclose(mod.eval_magnification_bias(*profile_pars, alpha=alpha),
+                        cfg['numcosmo_profiles']['mu']**(alpha-1), 1.e3*reltol)
+
         beta_s_mean = 0.6
         beta_s_square_mean = 0.4
         source_redshift_inf = 1000.
         gammat_inf = mod.eval_tangential_shear(
             profile_pars[0], profile_pars[1], source_redshift_inf)
         kappa_inf = mod.eval_convergence(profile_pars[0], profile_pars[1], source_redshift_inf)
+
+        # Validate reduced tangential shear
         assert_allclose(
-            mod.eval_reduced_tangential_shear(*profile_pars, 'applegate14',
-                                              beta_s_mean, beta_s_square_mean),
-            beta_s_mean * gammat_inf/(1.0 - beta_s_square_mean / beta_s_mean * kappa_inf),
+            mod.eval_reduced_tangential_shear(*profile_pars[:2],
+                                              (beta_s_mean, beta_s_square_mean),
+                                              'beta', 'order1'),
+            beta_s_mean * gammat_inf/(1.0 - beta_s_mean * kappa_inf),
             1.0e-10)
         assert_allclose(
-            mod.eval_reduced_tangential_shear(*profile_pars, 'schrabback18',
-                                              beta_s_mean, beta_s_square_mean),
-            (1. + (beta_s_square_mean / (beta_s_mean * beta_s_mean) - 1.) * beta_s_mean *
-             kappa_inf) * (beta_s_mean * gammat_inf / (1. - beta_s_mean * kappa_inf)),
+            mod.eval_reduced_tangential_shear(*profile_pars[:2],
+                                              (beta_s_mean, beta_s_square_mean),
+                                              'beta', 'order2'),
+            (1.+(beta_s_square_mean/(beta_s_mean*beta_s_mean)-1.)*beta_s_mean*kappa_inf) \
+            *(beta_s_mean*gammat_inf/(1.-beta_s_mean*kappa_inf)),
             1.0e-10)
 
         # Validate magnification
-        assert_allclose(mod.eval_magnification(*profile_pars),
-                        1./((1-kappa)**2-abs(gammat)**2), 1.0e-10)
-        assert_allclose(mod.eval_magnification(*profile_pars),
-                        cfg['numcosmo_profiles']['mu'], 1.e2*reltol)
+        alpha = -1.78
+        assert_allclose(
+            mod.eval_magnification(*profile_pars[:2],
+                                   (beta_s_mean, beta_s_square_mean),
+                                   'beta', 'order1'),
+            1 + 2*beta_s_mean*kappa_inf,
+            1.0e-10)
+        assert_allclose(
+            mod.eval_magnification(*profile_pars[:2],
+                                   (beta_s_mean, beta_s_square_mean),
+                                   'beta', 'order2'),
+            1 + 2*beta_s_mean*kappa_inf\
+              + 3*beta_s_square_mean*kappa_inf**2 + beta_s_square_mean*gammat_inf**2,
+            1.0e-10)
 
         # Validate magnification bias
-        alpha = -1.78
-        assert_allclose(mod.eval_magnification_bias(*profile_pars, alpha=alpha),
-                        1./((1-kappa)**2-abs(gammat)**2)**(alpha-1), 1.0e-10)
-        assert_allclose(mod.eval_magnification_bias(*profile_pars, alpha=alpha),
-                        cfg['numcosmo_profiles']['mu']**(alpha-1), 1.e3*reltol)
+        assert_allclose(
+            mod.eval_magnification_bias(*profile_pars[:2],
+                                        (beta_s_mean, beta_s_square_mean),
+                                        alpha,
+                                        'beta', 'order1'),
+            1 + (alpha-1)*(2*beta_s_mean*kappa_inf),
+            1.0e-10)
+        assert_allclose(
+            mod.eval_magnification_bias(*profile_pars[:2],
+                                        (beta_s_mean, beta_s_square_mean),
+                                        alpha,
+                                        'beta', 'order2'),
+            1 + (alpha-1)*(2*beta_s_mean*kappa_inf)\
+              + (alpha-1)*(beta_s_square_mean*gammat_inf**2)\
+              +(2*alpha-1)*(alpha-1)*beta_s_square_mean*kappa_inf**2,
+            1.0e-10)
+
+
 
         # Check that shear, reduced shear and convergence return zero
         # and magnification and magnification_bias return one
diff --git a/tests/test_theory_parent.py b/tests/test_theory_parent.py
index edfe2c274..53631e0e1 100644
--- a/tests/test_theory_parent.py
+++ b/tests/test_theory_parent.py
@@ -4,7 +4,6 @@
 from numpy.testing import assert_raises, assert_allclose, assert_equal
 import clmm.theory as theo
 from clmm.theory.parent_class import CLMModeling
-from clmm.utils import compute_beta_s_square_mean, compute_beta_s_mean
 
 def test_unimplemented(modeling_data):
     """ Unit tests abstract class unimplemented methdods """
@@ -26,11 +25,16 @@ def test_unimplemented(modeling_data):
     assert_raises(NotImplementedError, mod.eval_mean_surface_density, [0.3], 0.3)
     assert_raises(NotImplementedError, mod.eval_excess_surface_density, [0.3], 0.3)
     assert_raises(NotImplementedError, mod.eval_tangential_shear, [0.3], 0.3, 0.5)
-    assert_raises(NotImplementedError, mod.eval_reduced_tangential_shear, [0.3], 0.3, 0.5)
-    assert_raises(NotImplementedError, mod.eval_reduced_tangential_shear, [0.3], 0.3, 0.5, 'applegate14', 0.6, 0.4)
-    assert_raises(NotImplementedError, mod.eval_reduced_tangential_shear, [0.3], 0.3, 0.5, 'schrabback18', 0.6, 0.4)
     assert_raises(NotImplementedError, mod.eval_convergence, [0.3], 0.3, 0.5)
+    assert_raises(NotImplementedError, mod.eval_reduced_tangential_shear, [0.3], 0.3, 0.5)
+    assert_raises(NotImplementedError, mod.eval_reduced_tangential_shear, [0.3], 0.3, (0.6, 0.4),
+                  'beta', 'order1')
     assert_raises(NotImplementedError, mod.eval_magnification, [0.3], 0.3, 0.5)
+    assert_raises(NotImplementedError, mod.eval_magnification, [0.3], 0.3, (0.6, 0.4),
+                  'beta', 'order1')
+    assert_raises(NotImplementedError, mod.eval_magnification_bias, [0.3], 0.3, 0.5, 3.)
+    assert_raises(NotImplementedError, mod.eval_magnification_bias, [0.3], 0.3, (0.6, 0.4), 3.,
+                  'beta', 'order1')
 
 
 def test_instantiate(modeling_data):
@@ -104,43 +108,6 @@ def test_instantiate(modeling_data):
     reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, np.repeat(z_src, len(r_proj)))
     assert_allclose(reduced_shear, shear/(1.0-convergence), rtol=1.0e-12)
 
-    beta_s_mean = 0.9
-    beta_s_square_mean = 0.6
-    source_redshift_inf = 1000. 
-    shear_inf = mod.eval_tangential_shear(r_proj, z_cl, source_redshift_inf)
-    convergence_inf = mod.eval_convergence(r_proj, z_cl, source_redshift_inf)
-
-    #Tests with pre-fixed beta values
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, z_src, 'applegate14', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, beta_s_mean * shear_inf/(1.0 - beta_s_square_mean / beta_s_mean * convergence_inf), rtol=1.0e-12)
-
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, np.repeat(z_src, len(r_proj)), 'applegate14', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, beta_s_mean * shear_inf/(1.0 - beta_s_square_mean / beta_s_mean * convergence_inf), rtol=1.0e-12)
-
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, z_src, 'schrabback18', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, (1. + (beta_s_square_mean / (beta_s_mean * beta_s_mean) - 1.) * beta_s_mean * convergence_inf) * (beta_s_mean * shear_inf / (1. - beta_s_mean * convergence_inf)), rtol=1.0e-12)
-
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, np.repeat(z_src, len(r_proj)), 'schrabback18', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, (1. + (beta_s_square_mean / (beta_s_mean * beta_s_mean) - 1.) * beta_s_mean * convergence_inf) * (beta_s_mean * shear_inf / (1. - beta_s_mean * convergence_inf)), rtol=1.0e-12)
-
-    #Tests where the function computes the beta values
-    beta_s_mean = None
-    beta_s_square_mean = None
-
-    beta_s_square_test = compute_beta_s_square_mean(z_cl, 1000., mod.cosmo)
-    beta_s_test = compute_beta_s_mean(z_cl, 1000., mod.cosmo)   
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, z_src, 'applegate14', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, beta_s_test* shear_inf/(1.0 - beta_s_square_test / beta_s_test * convergence_inf), rtol=1.0e-12)
-
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, np.repeat(z_src, len(r_proj)), 'applegate14', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, beta_s_test * shear_inf/(1.0 - beta_s_square_test / beta_s_test * convergence_inf), rtol=1.0e-12)
-
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, z_src, 'schrabback18', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, (1. + (beta_s_square_test / (beta_s_test * beta_s_test) - 1.) * beta_s_test * convergence_inf) * (beta_s_test * shear_inf / (1. - beta_s_test * convergence_inf)), rtol=1.0e-12)
-
-    reduced_shear = mod.eval_reduced_tangential_shear(r_proj, z_cl, np.repeat(z_src, len(r_proj)), 'schrabback18', beta_s_mean, beta_s_square_mean)
-    assert_allclose(reduced_shear, (1. + (beta_s_square_test / (beta_s_test * beta_s_test) - 1.) * beta_s_test * convergence_inf) * (beta_s_test * shear_inf / (1. - beta_s_test * convergence_inf)), rtol=1.0e-12)
-
     assert_raises(ValueError, mod.eval_critical_surface_density, z_cl, use_pdz=True)
 
 def test_einasto(modeling_data):
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 5f0340773..2a17270ad 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -3,11 +3,12 @@
 import numpy as np
 from numpy.testing import assert_raises, assert_allclose, assert_equal
 from scipy.integrate import quad
+import clmm
 import clmm.utils as utils
-import clmm.theory as md
 from clmm.utils import (
     compute_radial_averages, make_bins, convert_shapes_to_epsilon, arguments_consistency,
     validate_argument)
+from clmm import z_distributions as zdist
 
 
 TOLERANCE = {'rtol': 1.0e-6, 'atol': 0}
@@ -223,7 +224,7 @@ def test_convert_units():
     """
     # Make an astropy cosmology object for testing
     # cosmo = FlatLambdaCDM(H0=70., Om0=0.3)
-    cosmo = md.Cosmology(H0=70.0, Omega_dm0=0.3-0.045, Omega_b0=0.045)
+    cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.3-0.045, Omega_b0=0.045)
 
     # Test that each unit is supported
     utils.convert_units(1.0, 'radians', 'degrees')
@@ -429,37 +430,45 @@ def test_validate_argument():
 
 def test_beta_functions():
     z_cl = 1.0
-    z_s = 2.4 
+    z_s = 2.4
     z_inf =1000.
     zmax = 15.0
     nsteps = 1000
     zmin = z_cl + 0.1
-    z_int = np.linspace(zmin, zmax, nsteps)    
-    cosmo = md.Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045,
+    z_int = np.linspace(zmin, zmax, nsteps)
+    cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045,
                   Omega_b0=0.045, Omega_k0=0.0)
-    beta_test = np.heaviside(z_s-z_cl, 0) * cosmo.eval_da_z1z2(z_cl, z_s) / cosmo.eval_da(z_s) 
+    beta_test = np.heaviside(z_s-z_cl, 0) * cosmo.eval_da_z1z2(z_cl, z_s) / cosmo.eval_da(z_s)
     beta_s_test = utils.compute_beta(z_s, z_cl, cosmo) / utils.compute_beta(z_inf, z_cl, cosmo)
 
-    def pdz(z):
-        return (z**1.24)*np.exp(-(z/0.51)**1.01)
-    
-    def integrand1(z_i, z_cl=z_cl, cosmo=cosmo):
-        return utils.compute_beta(z_i, z_cl, cosmo) * pdz(z_i)
-    
-    def integrand2(z_i, z_inf=z_inf, z_cl=z_cl, cosmo=cosmo):
-        return utils.compute_beta_s(z_i, z_cl, z_inf, cosmo) * pdz(z_i)
-    
-    def integrand3(z_i, z_inf=z_inf, z_cl=z_cl, cosmo=cosmo):
-        return utils.compute_beta_s(z_i, z_cl, z_inf, cosmo)**2 * pdz(z_i)
-
-    test1 = utils.compute_beta(z_s, z_cl, cosmo)
-    test2 = utils.compute_beta_s(z_s, z_cl, z_inf, cosmo)
-    test3 = utils.compute_beta_mean(z_cl, cosmo, zmax)
-    test4 = utils.compute_beta_s_mean(z_cl, z_inf,cosmo, zmax)
-    test5 = utils.compute_beta_s_square_mean(z_cl, z_inf,cosmo, zmax)
-    
-    assert_allclose(test1, beta_test, **TOLERANCE)
-    assert_allclose(test2, beta_s_test, **TOLERANCE)
-    assert_allclose(test3, quad(integrand1, zmin, zmax)[0] / quad(pdz, zmin, zmax)[0], **TOLERANCE)
-    assert_allclose(test4, quad(integrand2, zmin, zmax)[0] / quad(pdz, zmin, zmax)[0], **TOLERANCE)
-    assert_allclose(test5, quad(integrand3, zmin, zmax)[0] / quad(pdz, zmin, zmax)[0], **TOLERANCE)
+    assert_allclose(utils.compute_beta(z_s, z_cl, cosmo), beta_test, **TOLERANCE)
+    assert_allclose(utils.compute_beta_s(z_s, z_cl, z_inf, cosmo), beta_s_test, **TOLERANCE)
+
+
+    for model in (None, zdist.chang2013, zdist.desc_srd):
+
+        # None defaults to chang2013 for compute_beta* functions
+
+        test1 = utils.compute_beta_mean(z_cl, cosmo, zmax, z_distrib_func=model)
+        test2 = utils.compute_beta_s_mean(z_cl, z_inf, cosmo, zmax, z_distrib_func=model)
+        test3 = utils.compute_beta_s_square_mean(z_cl, z_inf, cosmo, zmax,
+                                                 z_distrib_func=model)
+
+        if model is None:
+            model = zdist.chang2013
+
+        def integrand1(z_i, z_cl=z_cl, cosmo=cosmo):
+            return utils.compute_beta(z_i, z_cl, cosmo) * model(z_i)
+
+        def integrand2(z_i, z_inf=z_inf, z_cl=z_cl, cosmo=cosmo):
+            return utils.compute_beta_s(z_i, z_cl, z_inf, cosmo) * model(z_i)
+
+        def integrand3(z_i, z_inf=z_inf, z_cl=z_cl, cosmo=cosmo):
+            return utils.compute_beta_s(z_i, z_cl, z_inf, cosmo)**2 * model(z_i)
+
+        assert_allclose(test1, quad(integrand1, zmin, zmax)[0] / quad(model, zmin, zmax)[0],
+                        **TOLERANCE)
+        assert_allclose(test2, quad(integrand2, zmin, zmax)[0] / quad(model, zmin, zmax)[0],
+                        **TOLERANCE)
+        assert_allclose(test3, quad(integrand3, zmin, zmax)[0] / quad(model, zmin, zmax)[0],
+                        **TOLERANCE)