diff --git a/clmm/__init__.py b/clmm/__init__.py
index adb901357..27802736c 100644
--- a/clmm/__init__.py
+++ b/clmm/__init__.py
@@ -12,6 +12,7 @@
     compute_convergence,
     compute_excess_surface_density,
     compute_excess_surface_density_2h,
+    compute_excess_surface_density_triaxial,
     compute_magnification,
     compute_magnification_bias,
     compute_magnification_bias_from_magnification,
@@ -27,4 +28,4 @@
 )
 from .utils import compute_radial_averages, convert_units, make_bins
 
-__version__ = "1.17.0"
+__version__ = "1.18.0"
diff --git a/clmm/clusterensemble.py b/clmm/clusterensemble.py
index d37591545..24d23e67f 100644
--- a/clmm/clusterensemble.py
+++ b/clmm/clusterensemble.py
@@ -31,14 +31,20 @@ class ClusterEnsemble:
 
         * "tan_sc" : tangential component computed with sample covariance
         * "cross_sc" : cross component computed with sample covariance
+        * "quad_4theta_sc" : quadrupole 4theta component computed with sample covariance
+        * "quad_const_sc" : quadrupole constant component computed with sample covariance
         * "tan_bs" : tangential component computed with bootstrap
         * "cross_bs" : cross component computed with bootstrap
+        * "quad_4theta_bs" : quadrupole 4theta component computed with bootstrap
+        * "quad_const_bs" : quadrupole constant component computed with bootstrap
         * "tan_jk" : tangential component computed with jackknife
         * "cross_jk" : cross component computed with jackknife
+        * "quad_4theta_jk" : quadrupole 4theta component computed with jackknife
+        * "quad_const_jk" : quadrupole constant component computed with jackknife
 
     """
 
-    def __init__(self, unique_id, gc_list=None, **kwargs):
+    def __init__(self, unique_id, gc_list=None, include_quadrupole=False, **kwargs):
         """Initializes a ClusterEnsemble object
         Parameters
         ----------
@@ -51,6 +57,7 @@ def __init__(self, unique_id, gc_list=None, **kwargs):
             raise TypeError(f"unique_id incorrect type: {type(unique_id)}")
         self.unique_id = unique_id
         self.data = GCData(meta={"bin_units": None, "radius_min": None, "radius_max": None})
+        self.include_quadrupole = include_quadrupole
         if gc_list is not None:
             self._add_values(gc_list, **kwargs)
             weights_out = kwargs.get("weights_out", "W_l")
@@ -63,10 +70,16 @@ def __init__(self, unique_id, gc_list=None, **kwargs):
         self.cov = {
             "tan_sc": None,
             "cross_sc": None,
+            "quad_4theta_sc": None,
+            "quad_const_sc": None,
             "tan_bs": None,
             "cross_bs": None,
+            "quad_4theta_bs": None,
+            "quad_const_bs": None,
             "tan_jk": None,
             "cross_jk": None,
+            "quad_4theta_jk": None,
+            "quad_const_jk": None,
         }
 
     def _add_values(self, gc_list, **kwargs):
@@ -90,6 +103,7 @@ def __len__(self):
         """Returns length of ClusterEnsemble"""
         return len(self.data)
 
+    # pylint: disable=R0914
     def make_individual_radial_profile(
         self,
         galaxycluster,
@@ -99,10 +113,16 @@ def make_individual_radial_profile(
         cosmo=None,
         tan_component_in="et",
         cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
         tan_component_out="gt",
         cross_component_out="gx",
+        quad_4theta_component_out="g_quad_4theta",
+        quad_const_component_out="g_quad_const",
         tan_component_in_err=None,
         cross_component_in_err=None,
+        quad_4theta_component_in_err=None,
+        quad_const_component_in_err=None,
         use_weights=True,
         weights_in="w_ls",
         weights_out="W_l",
@@ -137,18 +157,38 @@ def make_individual_radial_profile(
         cross_component_in: string, optional
             Name of the cross component column in `galcat` to be binned.
             Default: 'ex'
+        quad_4theta_component_in: string, optional
+            Name of the quadrupole 4theta component column in `galcat` to be binned.
+            Default: 'e_quad_4theta'
+        quad_const_component_in: string, optional
+            Name of the quadrupole constant component column in `galcat` to be binned.
+            Default: 'e_quad_const'
         tan_component_out: string, optional
             Name of the tangetial component binned column to be added in profile table.
             Default: 'gt'
         cross_component_out: string, optional
             Name of the cross component binned profile column to be added in profile table.
             Default: 'gx'
+        quad_4theta_component_out: string, optional
+            Name of the quadrupole 4theta component binned profile column to be added
+            in profile table.
+            Default: 'g_quad_4theta'
+        quad_const_component_out: string, optional
+            Name of the quadrupole constant component binned profile column to be added
+            in profile table.
+            Default: 'g_quad_const'
         tan_component_in_err: string, None, optional
             Name of the tangential component error column in `galcat` to be binned.
             Default: None
         cross_component_in_err: string, None, optional
             Name of the cross component error column in `galcat` to be binned.
             Default: None
+        quad_4theta_component_in_err: string, None, optional
+            Name of the quadrupole 4theta component error column in `galcat` to be binned.
+            Default: None
+        quad_const_component_in_err: string, None, optional
+            Name of the quadrupole constant component error column in `galcat` to be binned.
+            Default: None
         weights_in : str, None
             Name of the weight column in `galcat` to be considered in binning.
         weights_out : str, None
@@ -170,13 +210,14 @@ def make_individual_radial_profile(
         profile_table = galaxycluster.make_radial_profile(
             include_empty_bins=True, gal_ids_in_bins=False, add=False, **tb_kwargs
         )
-
         self.add_individual_radial_profile(
             galaxycluster,
             profile_table,
-            tan_component_out,
-            cross_component_out,
-            weights_out,
+            tan_component=tan_component_out,
+            cross_component=cross_component_out,
+            quad_4theta_component=quad_4theta_component_out,
+            quad_const_component=quad_const_component_out,
+            weights=weights_out,
         )
 
     def add_individual_radial_profile(
@@ -185,6 +226,8 @@ def add_individual_radial_profile(
         profile_table,
         tan_component="gt",
         cross_component="gx",
+        quad_4theta_component="g_quad_4theta",
+        quad_const_component="g_quad_const",
         weights="W_l",
     ):
         """Compute the individual shear profile from a single GalaxyCluster object
@@ -203,6 +246,12 @@ def add_individual_radial_profile(
         cross_component: string, optional
             Name of the cross component binned profile column in the profile table.
             Default: 'gx'
+        quad_4theta_component: string, optional
+            Name of the quadrupole 4theta component binned profile column in the profile table.
+            Default: 'g_quad_4theta'
+        quad_const_component: string, optional
+            Name of the quadrupole constant component binned profile column in the profile table.
+            Default: 'g_quad_const'
         weights : str, None
             Name of the weight binned column in the profile table.
         """
@@ -215,7 +264,11 @@ def add_individual_radial_profile(
         cl_cosmo = profile_table.meta.get("cosmo", None)
         self.data.update_info_ext_valid("cosmo", self.data, cl_cosmo, overwrite=False)
 
-        tbcols = ("radius", tan_component, cross_component, weights)
+        prof_cols = [tan_component, cross_component]
+        if self.include_quadrupole:
+            prof_cols += [quad_4theta_component, quad_const_component]
+        tbcols = ("radius", *prof_cols, weights)
+
         data_to_save = [
             galaxycluster.unique_id,
             galaxycluster.ra,
@@ -238,7 +291,14 @@ def _check_empty_data(self):
                 "for each cluster in your catalog"
             )
 
-    def make_stacked_radial_profile(self, tan_component="gt", cross_component="gx", weights="W_l"):
+    def make_stacked_radial_profile(
+        self,
+        tan_component="gt",
+        cross_component="gx",
+        quad_4theta_component="g_quad_4theta",
+        quad_const_component="g_quad_const",
+        weights="W_l",
+    ):
         """Computes stacked profile and mean separation distances and add it internally
         to `stacked_data`.
 
@@ -250,15 +310,25 @@ def make_stacked_radial_profile(self, tan_component="gt", cross_component="gx",
         cross_component : string, optional
             Name of the cross component column in `data`.
             Default: 'gx'
+        quad_4theta_component : string, optional
+            Name of the quadrupole 4theta component column in `data`.
+            Default: 'g_quad_4theta'
+        quad_const_component : string, optional
+            Name of the quadrupole constant component column in `data`.
+            Default: 'g_quad_const'
         weights : str
             Name of the weights column in `data`.
         """
         self._check_empty_data()
 
+        prof_cols = [tan_component, cross_component]
+        if self.include_quadrupole:
+            prof_cols += [quad_4theta_component, quad_const_component]
+
         radius, components = make_stacked_radial_profile(
             self.data["radius"],
             self.data[weights],
-            [self.data[tan_component], self.data[cross_component]],
+            [self.data[_col] for _col in prof_cols],
         )
         self.stacked_data = GCData(
             [
@@ -268,10 +338,16 @@ def make_stacked_radial_profile(self, tan_component="gt", cross_component="gx",
                 *components,
             ],
             meta={k: v for k, v in self.data.meta.items() if k not in ("radius_min", "radius_max")},
-            names=("radius_min", "radius_max", "radius", tan_component, cross_component),
+            names=("radius_min", "radius_max", "radius", *prof_cols),
         )
 
-    def compute_sample_covariance(self, tan_component="gt", cross_component="gx"):
+    def compute_sample_covariance(
+        self,
+        tan_component="gt",
+        cross_component="gx",
+        quad_4theta_component="g_quad_4theta",
+        quad_const_component="g_quad_const",
+    ):
         """Compute Sample covariance matrix for cross and tangential and cross
         stacked profiles and updates .cov dict (`tan_sc`, `cross_sc`).
 
@@ -283,15 +359,33 @@ def compute_sample_covariance(self, tan_component="gt", cross_component="gx"):
         cross_component : string, optional
             Name of the cross component column in `data`.
             Default: 'gx'
+        quad_4theta_component : string, optional
+            Name of the quadrupole 4theta component column in `data`.
+            Default: 'g_quad_4theta'
+        quad_const_component : string, optional
+            Name of the quadrupole constant component column in `data`.
+            Default: 'g_quad_const'
         """
         self._check_empty_data()
 
+        prof_cols = [tan_component, cross_component]
+        if self.include_quadrupole:
+            prof_cols += [quad_4theta_component, quad_const_component]
+
         n_catalogs = len(self.data)
-        self.cov["tan_sc"] = np.cov(self.data[tan_component].T, bias=False) / n_catalogs
-        self.cov["cross_sc"] = np.cov(self.data[cross_component].T, bias=False) / n_catalogs
+        for name, col in zip(
+            ("tan_sc", "cross_sc", "quad_4theta_sc", "quad_const_sc"),
+            prof_cols,
+        ):
+            self.cov[name] = np.cov(self.data[col].T, bias=False) / n_catalogs
 
     def compute_bootstrap_covariance(
-        self, tan_component="gt", cross_component="gx", n_bootstrap=10
+        self,
+        tan_component="gt",
+        cross_component="gx",
+        quad_4theta_component="g_quad_4theta",
+        quad_const_component="g_quad_const",
+        n_bootstrap=10,
     ):
         """Compute the bootstrap covariance matrix, add boostrap covariance matrix for
         tangential and cross stacked profiles and updates .cov dict (`tan_jk`, `cross_bs`).
@@ -304,31 +398,44 @@ def compute_bootstrap_covariance(
         cross_component : string, optional
             Name of the cross component column in `data`.
             Default: 'gx'
+        quad_4theta_component : string, optional
+            Name of the quadrupole 4theta component column in `data`.
+            Default: 'g_quad_4theta'
+        quad_const_component : string, optional
+            Name of the quadrupole constant component column in `data`.
+            Default: 'g_quad_const'
         n_bootstrap : int
             number of bootstrap resamplings
         """
         self._check_empty_data()
         n_catalogs = len(self)
+        cluster_index_bootstrap = np.random.choice(np.arange(n_catalogs), (n_bootstrap, n_catalogs))
 
-        cluster_index = np.arange(n_catalogs)
-        cluster_index_bootstrap = [
-            np.random.choice(cluster_index, n_catalogs) for n_boot in range(n_bootstrap)
-        ]
+        prof_cols = [tan_component, cross_component]
+        if self.include_quadrupole:
+            prof_cols += [quad_4theta_component, quad_const_component]
 
-        gt_boot, gx_boot = make_stacked_radial_profile(
+        g_boot_components = make_stacked_radial_profile(
             self["radius"][None, cluster_index_bootstrap][0].transpose(1, 2, 0),
             self["W_l"][None, cluster_index_bootstrap][0].transpose(1, 2, 0),
-            [
-                self[tan_component][None, cluster_index_bootstrap][0].transpose(1, 2, 0),
-                self[cross_component][None, cluster_index_bootstrap][0].transpose(1, 2, 0),
-            ],
+            [self[col][None, cluster_index_bootstrap][0].transpose(1, 2, 0) for col in prof_cols],
         )[1]
 
         coeff = (n_catalogs / (n_catalogs - 1)) ** 2
-        self.cov["tan_bs"] = coeff * np.cov(np.array(gt_boot), bias=False, ddof=0)
-        self.cov["cross_bs"] = coeff * np.cov(np.array(gx_boot), bias=False)
+        for name, component in zip(
+            ("tan_bs", "cross_bs", "quad_4theta_bs", "quad_const_bs"),
+            g_boot_components,
+        ):
+            self.cov[name] = coeff * np.cov(np.array(component), bias=False, ddof=0)
 
-    def compute_jackknife_covariance(self, tan_component="gt", cross_component="gx", n_side=16):
+    def compute_jackknife_covariance(
+        self,
+        tan_component="gt",
+        cross_component="gx",
+        quad_4theta_component="g_quad_4theta",
+        quad_const_component="g_quad_const",
+        n_side=16,
+    ):
         """Compute the jackknife covariance matrix, add boostrap covariance matrix for
         tangential and cross stacked profiles and updates .cov dict (`tan_jk`, `cross_jk`).
 
@@ -342,6 +449,12 @@ def compute_jackknife_covariance(self, tan_component="gt", cross_component="gx",
         cross_component : string, optional
             Name of the cross component column in `data`.
             Default: 'gx'
+        quad_4theta_component : string, optional
+            Name of the quadrupole 4theta component column in `data`.
+            Default: 'g_quad_4theta'
+        quad_const_component : string, optional
+            Name of the quadrupole constant component column in `data`.
+            Default: 'g_quad_const'
         n_side : int
             healpix sky area division parameter (number of sky area : 12*n_side^2)
         """
@@ -351,20 +464,30 @@ def compute_jackknife_covariance(self, tan_component="gt", cross_component="gx",
 
         pixels = healpy.ang2pix(n_side, self.data["ra"], self.data["dec"], nest=True, lonlat=True)
         pixels_list_unique = np.unique(pixels)
-        gt_jack, gx_jack = [], []
+
+        prof_cols = [tan_component, cross_component]
+        if self.include_quadrupole:
+            prof_cols += [quad_4theta_component, quad_const_component]
+
+        g_jack_components = [[] for _ in prof_cols]
         for hp_list_delete in pixels_list_unique:
             mask = ~np.isin(pixels, hp_list_delete)
-            gt, gx = make_stacked_radial_profile(
-                self["radius"][mask],
-                self["W_l"][mask],
-                [self[tan_component][mask], self[cross_component][mask]],
-            )[1]
-            gt_jack.append(gt)
-            gx_jack.append(gx)
+            for i, component in enumerate(
+                make_stacked_radial_profile(
+                    self["radius"][mask],
+                    self["W_l"][mask],
+                    [self[col][mask] for col in prof_cols],
+                )[1]
+            ):
+                g_jack_components[i].append(component)
+
         n_jack = pixels_list_unique.size
         coeff = (n_jack - 1) ** 2 / (n_jack)
-        self.cov["tan_jk"] = coeff * np.cov(np.transpose(gt_jack), bias=False, ddof=0)
-        self.cov["cross_jk"] = coeff * np.cov(np.transpose(gx_jack), bias=False, ddof=0)
+        for name, component in zip(
+            ("tan_jk", "cross_jk", "quad_4theta_jk", "quad_const_jk"),
+            g_jack_components,
+        ):
+            self.cov[name] = coeff * np.cov(np.transpose(component), bias=False, ddof=0)
 
     def save(self, filename, **kwargs):
         """Saves GalaxyCluster object to filename using Pickle"""
diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index 6af06c4c2..0cf8504fc 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -32,6 +32,8 @@ def compute_tangential_and_cross_components(
     geometry="curve",
     is_deltasigma=False,
     sigma_c=None,
+    include_quadrupole=False,
+    phi_major=None,
     validate_input=True,
 ):
     r"""Computes tangential- and cross- components for shear or ellipticity
@@ -63,14 +65,24 @@ def compute_tangential_and_cross_components(
         \tan\phi = & \frac{\delta_s-\delta_l}{\left(\alpha_l-\alpha_s\right)\cos(\delta_l)}
 
     The tangential, :math:`g_t`, and cross, :math:`g_x`, ellipticity/shear components are
-    calculated using the two ellipticity/shear components :math:`g_1` and :math:`g_2` of the
-    source galaxies, following Eq.7 and Eq.8 in Schrabback et al. (2018), arXiv:1611:03866 which
-    is consistent with arXiv:0509252
+    calculated using the two ellipticity/shear components :math:`g_1` and :math:`g_2` of the source
+    galaxies, following Eq.7 and Eq.8 in `Schrabback et al. 2017
+    <https://doi.org/10.1093/mnras/stx2666>`_ which is consistent with `Schneier 2006
+    <https://doi.org/10.1007/978-3-540-30310-7_3>`_:
 
     .. math::
         g_t =& -\left( g_1\cos\left(2\phi\right)+g_2\sin\left(2\phi\right)\right)\\
         g_x =& g_1 \sin\left(2\phi\right)-g_2\cos\left(2\phi\right)
 
+    The quadrupole ellipticity/shear components, :math:`g_4theta` and :math:`g_const`,
+    are also calculated using the two ellipticity/shear components :math:`g_1` and :math:`g_2`
+    of the source galaxies, following Eq.31 and Eq.34 in `Shin et al. (2018)
+    <https://doi.org/10.1093/mnras/stx3366>`_.
+
+    .. math::
+        g_{4\theta} =& \left( g_1\cos\left(4\phi\right)+g_2\sin\left(4\phi\right)\right)\\
+        g_{\rm const} =& g_1
+
     Finally, if  the critical surface density (:math:`\Sigma_\text{crit}`) is provided, an estimate
     of the excess surface density :math:`\widehat{\Delta\Sigma}` is obtained from
 
@@ -96,9 +108,10 @@ def compute_tangential_and_cross_components(
     shear2: array
         The measured shear (or reduced shear or ellipticity) of the source galaxies
     coordinate_system: str, optional
-        Coordinate system of the ellipticity components. Must be either 'celestial' or
-        euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
-        If not set, defaults to 'euclidean'.
+        Coordinate system of the ellipticity components. Must be either `celestial` or
+        `euclidean`. See `Rowe et al. 2015 <https://doi.org/10.48550/arXiv.1407.7676>`_ section 5.1
+        for more details.
+        Default is `euclidean`.
     geometry: str, optional
         Sky geometry to compute angular separation.
         Options are curve (uses astropy) or flat.
@@ -108,6 +121,16 @@ def compute_tangential_and_cross_components(
     sigma_c : None, array_like
         Critical (effective) surface density in units of :math:`M_\odot\ Mpc^{-2}`.
         Used only when is_deltasigma=True.
+    include_quadrupole: bool
+        If `True`, the quadrupole shear components (:math:`g_{4\theta},\; g_{\rm const}`) are
+        calculated .
+    phi_major : float, optional
+        Direction of the major axis of the input cluster in the unit of radian.
+        only needed when `include_quadrupole` is `True`.
+        This quantity is always in Euclidean coordinates, for celestial coordinates only set
+        the coordinate_system=`celestial`.
+        ``CLMM`` also provides a function to compute ``phi_major`` from the cluster members,
+        check :class:`clmm.data.calculate_major_axis` for details.
     validate_input: bool
         Validade each input argument
 
@@ -119,8 +142,13 @@ def compute_tangential_and_cross_components(
         Tangential shear (or assimilated quantity) for each source galaxy
     cross_component: array_like
         Cross shear (or assimilated quantity) for each source galaxy
+    IF include_quadrupole:
+        4theta_component: array_like
+            4theta shear component (or assimilated quantity) for each source galaxy
+        const_component: array_like
+            constant shear component (or assimilated quantity) for each source galaxy
     """
-    # pylint: disable-msg=too-many-locals
+    # pylint: disable-msg=too-many-locals,too-many-branches
     # Note: we make these quantities to be np.array so that a name is not passed from astropy
     # columns
     if coordinate_system is None:
@@ -135,6 +163,8 @@ def compute_tangential_and_cross_components(
         validate_argument(locals(), "shear1", "float_array")
         validate_argument(locals(), "shear2", "float_array")
         validate_argument(locals(), "geometry", str)
+        validate_argument(locals(), "is_deltasigma", bool)
+        validate_argument(locals(), "include_quadrupole", bool)
         validate_argument(locals(), "sigma_c", "float_array", none_ok=True)
         _validate_coordinate_system(locals(), "coordinate_system")
         ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency(
@@ -143,6 +173,7 @@ def compute_tangential_and_cross_components(
             prefix="Tangential- and Cross- shape components sources",
         )
         _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c)
+        validate_argument(locals(), "phi_major", float, none_ok=True)
     elif np.iterable(ra_source):
         ra_source_, dec_source_, shear1_, shear2_ = (
             np.array(col) for col in [ra_source, dec_source, shear1, shear2]
@@ -166,15 +197,29 @@ def compute_tangential_and_cross_components(
     else:
         raise NotImplementedError(f"Sky geometry {geometry} is not currently supported")
     # Compute the tangential and cross shears
-    tangential_comp = _compute_tangential_shear(shear1_, shear2_, phi)
-    cross_comp = _compute_cross_shear(shear1_, shear2_, phi)
+    output_components = [
+        _compute_tangential_shear(shear1_, shear2_, phi),
+        _compute_cross_shear(shear1_, shear2_, phi),
+    ]
 
+    if include_quadrupole:
+        if phi_major is None:
+            raise ValueError("phi_major should be provided for quadrupole computation")
+        if coordinate_system == "celestial":
+            phi_major = np.pi - phi_major
+        rotated_shear1, rotated_shear2 = _rotate_shear(shear1_, shear2_, phi_major)
+        # Compute the quadrupole shear components
+        output_components += [
+            _compute_4theta_shear(rotated_shear1, rotated_shear2, phi - phi_major),
+            rotated_shear1,  # const_comp
+        ]
+
+    # If the is_deltasigma flag is True, multiply the results by Sigma_crit.
     if sigma_c is not None:
         _sigma_c_arr = np.array(sigma_c)
-        tangential_comp *= _sigma_c_arr
-        cross_comp *= _sigma_c_arr
+        output_components = [component * _sigma_c_arr for component in output_components]
 
-    return angsep, tangential_comp, cross_comp
+    return angsep, *output_components
 
 
 def compute_background_probability(
@@ -369,9 +414,10 @@ def _compute_lensing_angles_flatsky(
     dec_source_list: array
         Declinations of each source galaxy in degrees
     coordinate_system: str, optional
-        Coordinate system of the ellipticity components. Must be either 'celestial' or
-        euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
-        Default is 'euclidean'.
+        Coordinate system of the ellipticity components. Must be either `celestial` or
+        `euclidean`. See `Rowe et al. 2015 <https://doi.org/10.48550/arXiv.1407.7676>`_ section 5.1
+        for more details.
+        Default is `euclidean`.
 
     Returns
     -------
@@ -419,9 +465,10 @@ def _compute_lensing_angles_astropy(
     dec_source_list: array
         Declinations of each source galaxy in degrees
     coordinate_system: str, optional
-        Coordinate system of the ellipticity components. Must be either 'celestial' or
-        euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
-        Default is 'euclidean'.
+        Coordinate system of the ellipticity components. Must be either `celestial` or
+        `euclidean`. See `Rowe et al. 2015 <https://doi.org/10.48550/arXiv.1407.7676>`_ section 5.1
+        for more details.
+        Default is `euclidean`.
 
     Returns
     -------
@@ -446,6 +493,76 @@ def _compute_lensing_angles_astropy(
     return angsep, phi
 
 
+def calculate_major_axis(ra_lens, dec_lens, ra_mem, dec_mem, weight_mem, remove_cg=True):
+    r"""Compute the major axis of a given cluster from the distribution of
+    its member galaxies using the position second moments.
+
+    The computation is done according to Eq. 5-10 of Shin et al. 2018, arXiv:1705.11167
+
+    Current implementation assumes the +RA direction is aligned with +x direction.
+
+    For extended descriptions of parameters, see `compute_shear()` documentation.
+
+    Parameters
+    ----------
+    ra_lens: float
+        Right ascension of the lensing cluster in degrees
+    dec_lens: float
+        Declination of the lensing cluster in degrees
+    ra_mem: array
+        Right ascensions of each cluster member galaxy in degrees
+    dec_mem: array
+        Declinations of each cluster member galaxy in degrees
+    weight_mem: array
+        Weight to be used for each cluster member galaxy
+    remove_cg: bool
+        Remove central galaxy from axis computation
+
+    Returns
+    -------
+    float
+        Major axis of cluster
+    """
+    sk_lens = SkyCoord(ra_lens * u.deg, dec_lens * u.deg, frame="icrs")
+    sk_mem = SkyCoord(np.array(ra_mem) * u.deg, np.array(dec_mem) * u.deg, frame="icrs")
+    # remove central galaxy
+    if remove_cg:
+        sk_mem = sk_mem[(np.array(ra_mem) != ra_lens) + (np.array(dec_mem) != dec_lens)]
+
+    # compute angles and weights
+    position_angle_mem = np.array(sk_lens.position_angle(sk_mem).radian + np.pi / 2.0)
+    separation_mem = np.array(sk_lens.separation(sk_mem).degree)
+
+    x_mem = separation_mem * np.cos(position_angle_mem)
+    y_mem = separation_mem * np.sin(position_angle_mem)
+
+    distance_weight_mem = 1.0 / (separation_mem**2)
+    weight_total_mem = np.array(weight_mem) * distance_weight_mem
+    weight_total_mem_norm = weight_total_mem / weight_total_mem.sum()
+
+    # Calcualte second moments of the member galaxies
+    ixx = (x_mem**2 * weight_total_mem_norm).sum()
+    iyy = (y_mem**2 * weight_total_mem_norm).sum()
+    ixy = (x_mem * y_mem * weight_total_mem_norm).sum()
+
+    # Transformation of the second moments to the direction of the major axis
+    return np.arctan2(2.0 * ixy, ixx - iyy) / 2.0
+
+
+def _rotate_shear(shear1, shear2, phi_major):
+    r"""Rotate shear components into the coordinate where +x axis is aligned with
+    the cluster major axis.
+
+    .. math::
+        g_rotated = g * \exp\left(-2i\phi\right)
+
+    For extended descriptions of parameters, see `compute_shear()` documentation.
+    """
+    shear_vec = shear1 + shear2 * 1.0j
+    rotated_shear_vec = shear_vec * np.exp(-2.0j * phi_major)
+    return np.real(rotated_shear_vec), np.imag(rotated_shear_vec)
+
+
 def _compute_tangential_shear(shear1, shear2, phi):
     r"""Compute the tangential shear given the two shears and azimuthal positions for
     a single source or list of sources.
@@ -475,6 +592,25 @@ def _compute_cross_shear(shear1, shear2, phi):
     return shear1 * np.sin(2.0 * phi) - shear2 * np.cos(2.0 * phi)
 
 
+def _compute_4theta_shear(shear1, shear2, phi):
+    r"""Compute the 4theta component of the quadrupole shear given the two shear components and
+    azimuthal positions for a single source or list of sources.
+
+    We compute the tangential shear following Eq. 31 of Shin et al. 2018, arXiv:1705.11167
+
+    .. math::
+        g_4theta = g_1\cos\left(4\phi\right)+g_2\sin\left(4\phi\right)
+
+    Note that here `\phi` is the position angle of the source galaxies
+    with respect to the major axis of the cluster.
+    Also, `shear1` and `shear2` are measured in the coordinate where
+    the +x axis is aligned with the major axis of the cluster.
+
+    For extended descriptions of parameters, see `compute_shear()' documentation.
+    """
+    return shear1 * np.cos(4.0 * phi) + shear2 * np.sin(4.0 * phi)
+
+
 def make_radial_profile(
     components,
     angsep,
@@ -545,9 +681,10 @@ def make_radial_profile(
         Array of individual galaxy weights. If specified, the radial binned profile is
         computed using a weighted average
     coordinate_system: str, optional
-        Coordinate system of the ellipticity components. Must be either 'celestial' or
-        euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
-        Default is 'euclidean'.
+        Coordinate system of the ellipticity components. Must be either `celestial` or
+        `euclidean`. See `Rowe et al. 2015 <https://doi.org/10.48550/arXiv.1407.7676>`_ section 5.1
+        for more details.
+        Default is `euclidean`.
     empty_bins_value: float, None
         Values to be assigned to empty bins.
 
@@ -656,7 +793,7 @@ def make_stacked_radial_profile(angsep, weights, components):
     Parameters
     ----------
     angsep: 2d array
-        Transvesal distances corresponding to each object with shape `n_obj, n_rad_bins`.
+        Transversal distances corresponding to each object with shape `n_obj, n_rad_bins`.
     weights: 2d array
         Weights corresponding to each objects with shape `n_obj, n_rad_bins`.
     components: list of 2d arrays
diff --git a/clmm/galaxycluster.py b/clmm/galaxycluster.py
index 4cb52d00a..4ba7a9040 100644
--- a/clmm/galaxycluster.py
+++ b/clmm/galaxycluster.py
@@ -6,6 +6,7 @@
 import warnings
 
 from .dataops import (
+    calculate_major_axis,
     compute_background_probability,
     compute_galaxy_weights,
     compute_tangential_and_cross_components,
@@ -35,19 +36,28 @@ class GalaxyCluster:
         Declination of galaxy cluster center (in degrees)
     z : float
         Redshift of galaxy cluster center
+    phi_major : float, None
+        Direction of the major axis of the cluster in radians, only used for triaxiality.
     galcat : GCData
         Table of background galaxy data containing at least galaxy_id, ra, dec, e1, e2, z
     validate_input: bool
         Validade each input argument
+    include_quadrupole: bool
+        If quadrupole WL be calculated.
     """
 
-    def __init__(self, *args, validate_input=True, **kwargs):
+    # pylint: disable=too-many-instance-attributes
+    # Eight is reasonable in this case.
+
+    def __init__(self, *args, validate_input=True, include_quadrupole=False, **kwargs):
         self.unique_id = None
         self.ra = None
         self.dec = None
         self.z = None
         self.galcat = None
+        self.phi_major = None
         self.validate_input = validate_input
+        self.include_quadrupole = include_quadrupole
         if len(args) > 0 or len(kwargs) > 0:
             self._add_values(*args, **kwargs)
             self._check_types()
@@ -60,6 +70,7 @@ def _add_values(
         dec: float,
         z: float,
         galcat: GCData,
+        phi_major: float = None,
     ):
         """Add values for all attributes"""
         self.unique_id = unique_id
@@ -67,6 +78,7 @@ def _add_values(
         self.dec = dec
         self.z = z
         self.galcat = galcat
+        self.phi_major = phi_major
 
     def _check_types(self):
         """Check types of all attributes"""
@@ -123,6 +135,25 @@ def _repr_html_(self):
             f"<br>{self.galcat._html_table()}"
         )
 
+    def set_phi_major(self, phi_major=None, info_mem=None):
+        r"""Sets the phi_major value for the cluster, is
+        obtained either from ``phi_major`` or from ``info_mem``.
+
+        Parameters
+        ----------
+        phi_major: float, optional
+            Cluster major axis direction (in radian with respect to +x).
+        info_mem: tuple, optional
+            (RA, DEC, weights) of cluster member galaxies as a tuple of array.
+        """
+        if phi_major is None:
+            if info_mem is None:
+                raise TypeError("Either phi_major or info_mem should be provided.")
+            phi_major = calculate_major_axis(self.ra, self.dec, *info_mem)
+        elif info_mem is not None:
+            raise TypeError("Only phi_major OR info_mem should be provided.")
+        self.phi_major = phi_major
+
     def add_critical_surface_density(self, cosmo, use_pdz=False, force=False):
         r"""Computes the critical surface density for each galaxy in `galcat`.
         It only runs if input cosmo != galcat cosmo or if `sigma_c` not in `galcat`.
@@ -217,6 +248,8 @@ def compute_tangential_and_cross_components(
         shape_component2="e2",
         tan_component="et",
         cross_component="ex",
+        quad_4theta_component="e_quad_4theta",
+        quad_const_component="e_quad_const",
         geometry="curve",
         is_deltasigma=False,
         use_pdz=False,
@@ -235,6 +268,8 @@ def compute_tangential_and_cross_components(
         geometry: `input` geometry
         is_deltasigma: `input` is_deltasigma
         coordinate_system: `galcat` coordinate_system
+        include_quadrupole: `input` include_quadrupole
+        phi_major: `cluster` major axis direction (in radian with respect to +x)
 
         Parameters
         ----------
@@ -252,6 +287,16 @@ def compute_tangential_and_cross_components(
             Name of the column to be added to the `galcat` astropy table that will contain the
             cross component computed from columns `shape_component1` and `shape_component2`.
             Default: `ex`
+        quad_4theta_component: string, optional
+            Name of the column to be added to the `galcat` astropy table that will contain the
+            4theta quarupole component computed from columns `shape_component1`
+            and `shape_component2`.
+            Default: `e_quad_4theta`
+        quad_const_component: string, optional
+            Name of the column to be added to the `galcat` astropy table that will contain the
+            constant quarupole component computed from columns `shape_component1`
+            and `shape_component2`.
+            Default: `e_quad_const`
         geometry: str, optional
             Sky geometry to compute angular separation.
             Options are curve (uses astropy) or flat.
@@ -259,7 +304,7 @@ def compute_tangential_and_cross_components(
             If `True`, the tangential and cross components returned are multiplied by Sigma_crit.
             Results in units of :math:`M_\odot\ Mpc^{-2}`
         cosmo: astropy cosmology object
-            Specifying a cosmology is required if `is_deltasigma` is True
+            Specifying a cosmology is required if `is_deltasigma` is `True`
         add: bool
             If `True`, adds the computed shears to the `galcat`
 
@@ -271,7 +316,17 @@ def compute_tangential_and_cross_components(
             Tangential shear (or assimilated quantity) for each source galaxy
         cross_component: array_like
             Cross shear (or assimilated quantity) for each source galaxy
+        quad_4theta_component: array_like
+            4theta quadrupole shear (or assimilated quantity) for each source galaxy
+            Returned only if include_quadrupole is `True`.
+        quad_const_component: array_like
+            constatnt quadrupole shear (or assimilated quantity) for each source galaxy
+            Returned only if include_quadrupole is `True`.
         """
+        # Check quadrupole input
+        if self.include_quadrupole and self.phi_major is None:
+            raise ValueError("Cluster phi_major is not set, run set_phi_major before.")
+
         # Check is all the required data is available
         col_dict = {
             "ra_source": "ra",
@@ -285,24 +340,33 @@ def compute_tangential_and_cross_components(
         cols = self._get_input_galdata(col_dict)
 
         # compute shears
-        angsep, tangential_comp, cross_comp = compute_tangential_and_cross_components(
+        angsep_and_components = compute_tangential_and_cross_components(
             is_deltasigma=is_deltasigma,
             ra_lens=self.ra,
             dec_lens=self.dec,
             geometry=geometry,
             validate_input=self.validate_input,
+            include_quadrupole=self.include_quadrupole,
+            phi_major=self.phi_major,
             coordinate_system=self.galcat.meta["coordinate_system"],
             **cols,
         )
+
         if add:
-            self.galcat["theta"] = angsep
-            self.galcat[tan_component] = tangential_comp
-            self.galcat[cross_component] = cross_comp
+
+            prof_cols = [tan_component, cross_component]
+            if self.include_quadrupole:
+                prof_cols += [quad_4theta_component, quad_const_component]
+
+            self.galcat["theta"] = angsep_and_components[0]
+            for _col, _comp_val in zip(prof_cols, angsep_and_components[1:]):
+                self.galcat[_col] = _comp_val
             if is_deltasigma:
                 sigmac_type = "effective" if use_pdz else "standard"
-                self.galcat.meta[f"{tan_component}_sigmac_type"] = sigmac_type
-                self.galcat.meta[f"{cross_component}_sigmac_type"] = sigmac_type
-        return angsep, tangential_comp, cross_comp
+                for _col in prof_cols:
+                    self.galcat.meta[f"{_col}_sigmac_type"] = sigmac_type
+
+        return angsep_and_components
 
     def compute_background_probability(
         self, use_pdz=False, add=True, p_background_name="p_background"
@@ -479,10 +543,16 @@ def make_radial_profile(
         cosmo=None,
         tan_component_in="et",
         cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
         tan_component_out="gt",
         cross_component_out="gx",
+        quad_4theta_component_out="g_quad_4theta",
+        quad_const_component_out="g_quad_const",
         tan_component_in_err=None,
         cross_component_in_err=None,
+        quad_4theta_component_in_err=None,
+        quad_const_component_in_err=None,
         include_empty_bins=False,
         gal_ids_in_bins=False,
         add=True,
@@ -498,7 +568,8 @@ def make_radial_profile(
         tangential shears or ellipticities and angular separation of the source galaxies
 
         Calls `clmm.dataops.make_radial_profile` with the following arguments:
-        components: `galcat` components (tan_component_in, cross_component_in, z)
+        components: `galcat` components (tan_component_in, cross_component_in, z) OR
+        (quad_4theta_component_in, quad_const_component_in, z) IF include_quadrupole
         angsep: `galcat` theta
         angsep_units: 'radians'
         bin_units: `input` bin_units
@@ -532,18 +603,37 @@ def make_radial_profile(
         cross_component_in: string, optional
             Name of the cross component column in `galcat` to be binned.
             Default: 'ex'
+        quad_4theta_component_in: string, optional
+            Name of the 4theta quadrupole component column in `galcat` to be binned.
+            Default: 'e_quad_4theta'
+        quad_const_component_in: string, optional
+            Name of the constant quadrupole component column in `galcat` to be binned.
+            Default: 'e_quad_const'
         tan_component_out: string, optional
             Name of the tangetial component binned column to be added in profile table.
             Default: 'gt'
         cross_component_out: string, optional
             Name of the cross component binned profile column to be added in profile table.
             Default: 'gx'
+        quad_4theta_component_out: string, optional
+            Name of the 4theta quadrupole component binned column to be added in profile table.
+            Default: 'g_quad_4theta'
+        quad_const_component_out: string, optional
+            Name of the constant quadrupole component binned profile column
+            to be added in profile table.
+            Default: 'g_quad_const'
         tan_component_in_err: string, None, optional
             Name of the tangential component error column in `galcat` to be binned.
             Default: None
         cross_component_in_err: string, None, optional
             Name of the cross component error column in `galcat` to be binned.
             Default: None
+        quad_4theta_component_in_err: string, None, optional
+            Name of the 4theta quadrupole component error column in `galcat` to be binned.
+            Default: None
+        quad_const_component_in_err: string, None, optional
+            Name of the constant quadrupole component error column in `galcat` to be binned.
+            Default: None
         include_empty_bins: bool, optional
             Also include empty bins in the returned table
         gal_ids_in_bins: bool, optional
@@ -572,10 +662,14 @@ def make_radial_profile(
         """
         # Too many local variables (19/15)
         # pylint: disable=R0914
-
-        if not all(
-            t_ in self.galcat.columns for t_ in (tan_component_in, cross_component_in, "theta")
-        ):
+        prof_col_in = [tan_component_in, cross_component_in]
+        prof_col_in_err = [tan_component_in_err, cross_component_in_err]
+        prof_col_out = [tan_component_out, cross_component_out]
+        if self.include_quadrupole:
+            prof_col_in += [quad_4theta_component_in, quad_const_component_in]
+            prof_col_in_err += [quad_4theta_component_in_err, quad_const_component_in_err]
+            prof_col_out += [quad_4theta_component_out, quad_const_component_out]
+        if not all(t_ in self.galcat.columns for t_ in (*prof_col_in, "theta")):
             raise TypeError(
                 "Shear or ellipticity information is missing. Galaxy catalog must have tangential"
                 "and cross shears (gt, gx) or ellipticities (et, ex). "
@@ -585,7 +679,7 @@ def make_radial_profile(
             raise TypeError("Missing galaxy redshifts!")
         # Compute the binned averages and associated errors
         profile_table, binnumber = make_radial_profile(
-            [self.galcat[n].data for n in (tan_component_in, cross_component_in, "z")],
+            [self.galcat[n].data for n in (*prof_col_in, "z")],
             angsep=self.galcat["theta"],
             angsep_units="radians",
             bin_units=bin_units,
@@ -597,14 +691,13 @@ def make_radial_profile(
             z_lens=self.z,
             validate_input=self.validate_input,
             components_error=[
-                None if n is None else self.galcat[n].data
-                for n in (tan_component_in_err, cross_component_in_err, None)
+                None if n is None else self.galcat[n].data for n in (*prof_col_in_err, None)
             ],
             weights=self.galcat[weights_in].data if use_weights else None,
             coordinate_system=self.galcat.meta["coordinate_system"],
         )
         # Reaname table columns
-        for i, name in enumerate([tan_component_out, cross_component_out, "z"]):
+        for i, name in enumerate([*prof_col_out, "z"]):
             profile_table.rename_column(f"p_{i}", name)
             profile_table.rename_column(f"p_{i}_err", f"{name}_err")
         # Reaname weights columns
@@ -638,6 +731,10 @@ def plot_profiles(
         tangential_component_error="gt_err",
         cross_component="gx",
         cross_component_error="gx_err",
+        quad_4theta_component="g_quad_4theta",
+        quad_4theta_component_error="g_quad_4theta_err",
+        quad_const_component="g_quad_const",
+        quad_const_component_error="g_quad_const_err",
         table_name="profile",
         xscale="linear",
         yscale="linear",
@@ -658,6 +755,19 @@ def plot_profiles(
         cross_component_error: str, optional
             Name of the column in the galcat Table corresponding to the uncertainty in the cross
             component of the shear or reduced shear. Default: 'gx_err'
+        quad_4theta_component: str, optional
+            Name of the column in the galcat Table corresponding to the 4theta quadrupole component
+            of the shear or reduced shear (Delta Sigma not yet implemented).
+            Default: 'g_quad_4theta'
+        quad_4theta_component_error: str, optional
+            Name of the column in the galcat Table corresponding to the uncertainty in
+            4theta quadrupole component of the shear or reduced shear. Default: 'g_quad_4theta_err'
+        quad_const_component: str, optional
+            Name of the column in the galcat Table corresponding to the constant quadrupole
+            component of the shear or reduced shear. Default: 'gconst'
+        quad_const_component_error: str, optional
+            Name of the column in the galcat Table corresponding to the uncertainty in the
+            constant quadrupole component of the shear or reduced shear. Default: 'gconst_er'
         table_name: str, optional
             Name of the GalaxyCluster() `.profile` attribute. Default: 'profile'
         xscale:
@@ -675,29 +785,38 @@ def plot_profiles(
         if not hasattr(self, table_name):
             raise ValueError(f"GalaxyClusters does not have a '{table_name}' table.")
         profile = getattr(self, table_name)
-        for col in (tangential_component, cross_component):
-            if col not in profile.columns:
-                raise ValueError(f"Column for plotting '{col}' does not exist.")
-        for col in (tangential_component_error, cross_component_error):
-            if col not in profile.columns:
-                warnings.warn(f"Column for plotting '{col}' does not exist.")
-        return plot_profiles(
-            rbins=profile["radius"],
-            r_units=profile.meta["bin_units"],
-            tangential_component=profile[tangential_component],
-            tangential_component_error=(
-                profile[tangential_component_error]
-                if tangential_component_error in profile.columns
-                else None
-            ),
-            cross_component=profile[cross_component],
-            cross_component_error=(
-                profile[cross_component_error] if cross_component_error in profile.columns else None
-            ),
-            xscale=xscale,
-            yscale=yscale,
-            tangential_component_label=tangential_component,
-            cross_component_label=cross_component,
+
+        prof_cols = [tangential_component, cross_component]
+        prof_cols_err = [tangential_component_error, cross_component_error]
+        if self.include_quadrupole:
+            prof_cols += [quad_4theta_component, quad_const_component]
+            prof_cols_err += [quad_4theta_component_error, quad_const_component_error]
+
+        for col in filter(lambda col: col not in profile.columns, prof_cols):
+            raise ValueError(f"Column for plotting '{col}' does not exist.")
+        for col in filter(lambda col: col not in profile.columns, prof_cols_err):
+            warnings.warn(f"Column for plotting '{col}' does not exist.")
+
+        return tuple(
+            plot_profiles(
+                rbins=profile["radius"],
+                r_units=profile.meta["bin_units"],
+                tangential_component=profile[comp1],
+                tangential_component_error=(
+                    profile[comp1_err] if comp1_err in profile.columns else None
+                ),
+                cross_component=profile[comp2],
+                cross_component_error=(
+                    profile[comp2_err] if comp2_err in profile.columns else None
+                ),
+                xscale=xscale,
+                yscale=yscale,
+                tangential_component_label=comp1,
+                cross_component_label=comp2,
+            )
+            for comp1, comp2, comp1_err, comp2_err in zip(
+                prof_cols[::2], prof_cols[1::2], prof_cols_err[::2], prof_cols_err[1::2]
+            )
         )
 
     def set_ra_lower(self, ra_low=0):
diff --git a/clmm/theory/__init__.py b/clmm/theory/__init__.py
index 10420f6fc..a037bd0d9 100644
--- a/clmm/theory/__init__.py
+++ b/clmm/theory/__init__.py
@@ -15,6 +15,7 @@
     compute_critical_surface_density_eff,
     compute_excess_surface_density,
     compute_excess_surface_density_2h,
+    compute_excess_surface_density_triaxial,
     compute_magnification,
     compute_magnification_bias,
     compute_mean_surface_density,
diff --git a/clmm/theory/func_layer.py b/clmm/theory/func_layer.py
index 9373d675c..f25ff8f99 100644
--- a/clmm/theory/func_layer.py
+++ b/clmm/theory/func_layer.py
@@ -1344,3 +1344,157 @@ def compute_magnification_bias(
 
     _modeling_object.validate_input = True
     return magnification_bias
+
+
+def compute_excess_surface_density_triaxial(
+    r_proj,
+    mdelta,
+    cdelta,
+    z_cl,
+    ell,
+    cosmo,
+    term,
+    n_grid=10000,
+    delta_mdef=200,
+    halo_profile_model="nfw",
+    massdef="mean",
+    alpha_ein=None,
+    r_mis=None,
+    mis_from_backend=False,
+    verbose=False,
+    use_projected_quad=False,
+    validate_input=True,
+):
+    r"""Compute the excess surface density lensing profile for the monopole, 4theta quadrupole,
+    or constant quadrupole component given in
+    `Shin et al. 2018 <https://doi.org/10.1093/mnras/stx3366>`_:
+
+    .. math::
+        \Delta\Sigma^{\rm mono} (R) \equiv \frac{2}{R^2}\int_0^R \mathrm{d} R' \Sigma_0(R') -
+        \Sigma_0(R)
+
+    .. math::
+        \Delta\Sigma^{4\theta} (R) \equiv \frac{\Sigma_{\rm crit}}{2\pi} \int \mathrm{d}\theta
+        \left( \gamma_1(R,\theta)\cos{4\theta} + \gamma_2(R,\theta) \sin{4\theta} \right) \\ =
+        \frac{\Sigma_2}{2} - \frac{3}{R^4}\int_0^R \mathrm{d}R' R'^3 \Sigma_2 (R')
+
+    .. math::
+        \Delta\Sigma^{\rm const} (R) \equiv \frac{\Sigma_{\rm crit}}{2\pi} \int \mathrm{d}\theta
+        \gamma_1(R,\theta) = \frac{\Sigma_2}{2} - \int_R^{\infty} \mathrm{d} R'
+        \frac{\Sigma_2(R')}{R'}
+
+    where :math:`\Sigma_0,\; \Sigma_2` come from multipole coefficients of :math:`\Sigma(R,
+    \theta)`:
+
+    .. math::
+        \Sigma_0(R) = \Sigma_{\rm sph}(R)\left[1 + \dfrac{\epsilon^2}{2} \left(\eta(R) +
+        \dfrac{\eta(R)^2}{2} + \frac{1}{2}\dfrac{\mathrm{d} \eta(R)}{\mathrm{d}\ln R} \right)
+        \right] + \mathcal{O}( \epsilon^3 )
+
+    .. math::
+        \Sigma_2(R) = - \epsilon \eta(R) \Sigma_{\rm sph} (R) \cos({2\theta_0}) + \mathcal{O}(
+        \epsilon^3 ),
+
+
+    with:
+
+    .. math::
+        \eta(R) = \frac{\mathrm{d} \ln{\Sigma_{\rm sph}(R)}} {\mathrm{d} \ln{R}},\;
+        \epsilon = \dfrac{1-q}{1+q}
+
+    and :math:`\theta_0` being the polar angle of the major axis (:math:`\theta_0=0` when we
+    align the major axis to the :math:`x`-axis).
+
+    Parameters
+    ----------
+    r_source : array
+        Radial distance of each source galaxy
+    mdelta : float
+        Mass of lens cluster
+    cdelta : array
+        Concentration of lens cluster
+    z_cluster : float
+        Redshift of lens cluster
+    ell : float
+        Ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
+        q=b/a (Ratio of major axis to the minor axis lengths)
+    cosmo : clmm.cosmology.Cosmology object
+        CLMM Cosmology object
+    term : str
+        The component of the Taylor expansion to return as described in Shin et al. 2018.
+        Must be in : 'mono', 'quad_4theta', 'quad_const'
+        - 'mono' : the ellipticity corrected monopole term
+        - 'quad_4theta' : the 4theta component of the quadrupole term
+        - 'quad_const' : the constant component of the quadrupole term
+    n_grid : int
+        Grid resolution for functions to be computed on.
+        Too low n_grid can lead to large deviations.
+    delta_mdef : int, optional
+        Mass overdensity definition; defaults to 200.
+    halo_profile_model : str, optional
+        Profile model parameterization (letter case independent):
+
+            * ``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``
+
+    alpha_ein : float, None, optional
+        If ``halo_profile_model=='einasto'``, set the value of the Einasto slope.
+        Option only available for the NumCosmo and CCL backends.
+        If None, use the default value of the backend. (0.25 for the NumCosmo backend and a
+        cosmology-dependent value for the CCL backend.)
+
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`
+    mis_from_backend : bool, optional
+        If True, use the projected surface density from the backend for miscentering
+        calculations. If False, use the (faster) CLMM exact analytical
+        implementation instead. (Default: False)
+    verbose : boolean, optional
+        If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL
+        backends.
+    use_projected_quad : bool
+        Only available for Einasto profile with CCL as the backend. If True, CCL will use
+        quad_vec instead of default FFTLog to calculate the surface density profile.
+        Default: False
+    validate_input : bool, optional
+        If True (default), the types of the arguments are checked before proceeding.
+
+    Returns
+    -------
+    dsmono : array
+        Component of delta sigma excess for the elliptical halo specified
+    """
+
+    _modeling_object.validate_input = validate_input
+    _modeling_object.set_cosmo(cosmo)
+    _modeling_object.set_halo_density_profile(
+        halo_profile_model=halo_profile_model, massdef=massdef, delta_mdef=delta_mdef
+    )
+    _modeling_object.set_concentration(cdelta)
+    _modeling_object.set_mass(mdelta)
+    if halo_profile_model == "einasto" or alpha_ein is not None:
+        _modeling_object.set_einasto_alpha(alpha_ein)
+    if halo_profile_model == "einasto" and _modeling_object.backend == "ccl":
+        _modeling_object.set_projected_quad(use_projected_quad)
+
+    delta_sigma = _modeling_object.eval_excess_surface_density_triaxial(
+        r_proj,
+        z_cl,
+        ell,
+        term,
+        r_mis=r_mis,
+        mis_from_backend=mis_from_backend,
+        verbose=verbose,
+        n_grid=n_grid,
+    )
+
+    _modeling_object.validate_input = True
+    return delta_sigma
diff --git a/clmm/theory/parent_class.py b/clmm/theory/parent_class.py
index f7b8e6097..6b78f2803 100644
--- a/clmm/theory/parent_class.py
+++ b/clmm/theory/parent_class.py
@@ -18,7 +18,7 @@
     compute_for_good_redshifts,
     validate_argument,
 )
-from . import miscentering
+from . import miscentering, triaxiality
 from .generic import (
     compute_magnification_bias_from_magnification,
     compute_profile_mass_in_radius,
@@ -278,6 +278,26 @@ def _eval_excess_surface_density_2h(
             2, r_proj, z_cl, halobias, logkbounds, ksteps, loglbounds, lsteps
         )
 
+    def _eval_excess_surface_density_triaxial(
+        self, surface_density_func, r_proj, z_cl, ell, term, n_grid=10000
+    ):
+        """eval individual terms of  excess surface density"""
+
+        args = (surface_density_func, r_proj, z_cl, ell, n_grid)
+
+        if term == "mono":
+            delta_sigma = self._eval_excess_surface_density(
+                r_proj, z_cl
+            ) + triaxiality.excess_surface_density_mono_correction(*args)
+        elif term == "quad_4theta":
+            delta_sigma = triaxiality.excess_surface_density_quad_4theta(*args)
+        elif term == "quad_const":
+            delta_sigma = triaxiality.excess_surface_density_quad_const(*args)
+        else:
+            raise ValueError(f"Unsupported term (='{term}')")
+
+        return delta_sigma
+
     def _eval_rdelta(self, z_cl):
         delta_mdef = self._get_delta_mdef(z_cl)
         return compute_rdelta(self.mdelta, z_cl, self.cosmo, self.massdef, delta_mdef)
@@ -684,8 +704,7 @@ def eval_surface_density(self, r_proj, z_cl, r_mis=None, mis_from_backend=False,
         if self.validate_input:
             validate_argument(locals(), "r_proj", "float_array", argmin=0)
             validate_argument(locals(), "z_cl", float, argmin=0)
-            if r_mis is not None:
-                validate_argument(locals(), "r_mis", float, argmin=0, eqmin=True)
+            validate_argument(locals(), "r_mis", float, argmin=0, eqmin=True, none_ok=True)
 
         if self.halo_profile_model == "einasto" and verbose:
             print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
@@ -725,8 +744,7 @@ def eval_mean_surface_density(
         if self.validate_input:
             validate_argument(locals(), "r_proj", "float_array", argmin=0)
             validate_argument(locals(), "z_cl", float, argmin=0)
-            if r_mis is not None:
-                validate_argument(locals(), "r_mis", float, argmin=0)
+            validate_argument(locals(), "r_mis", float, argmin=0, eqmin=True, none_ok=True)
 
         if self.halo_profile_model == "einasto" and verbose:
             print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
@@ -767,8 +785,7 @@ def eval_excess_surface_density(
         if self.validate_input:
             validate_argument(locals(), "r_proj", "float_array", argmin=0)
             validate_argument(locals(), "z_cl", float, argmin=0)
-            if r_mis is not None:
-                validate_argument(locals(), "r_mis", float, argmin=0, eqmin=True)
+            validate_argument(locals(), "r_mis", float, argmin=0, eqmin=True, none_ok=True)
 
         if self.halo_profile_model == "einasto" and verbose:
             print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
@@ -885,6 +902,112 @@ def eval_surface_density_2h(
             r_proj, z_cl, halobias, logkbounds, ksteps, loglbounds, lsteps
         )
 
+    def eval_excess_surface_density_triaxial(
+        self,
+        r_proj,
+        z_cl,
+        ell,
+        term,
+        r_mis=None,
+        mis_from_backend=False,
+        verbose=False,
+        n_grid=10000,
+    ):
+        r"""Compute the excess surface density lensing profile for the monopole, 4theta quadrupole,
+        or constant quadrupole component given in `Shin et al. 2018
+        <https://doi.org/10.1093/mnras/stx3366>`_:
+
+        .. math::
+            \Delta\Sigma^{\rm mono} (R) \equiv \frac{2}{R^2}\int_0^R \mathrm{d} R' \Sigma_0(R') -
+            \Sigma_0(R)
+
+        .. math::
+            \Delta\Sigma^{4\theta} (R) \equiv \frac{\Sigma_{\rm crit}}{2\pi} \int \mathrm{d}\theta
+            \left( \gamma_1(R,\theta)\cos{4\theta} + \gamma_2(R,\theta) \sin{4\theta} \right) \\ =
+            \frac{\Sigma_2}{2} - \frac{3}{R^4}\int_0^R \mathrm{d}R' R'^3 \Sigma_2 (R')
+
+        .. math::
+            \Delta\Sigma^{\rm const} (R) \equiv \frac{\Sigma_{\rm crit}}{2\pi} \int \mathrm{d}\theta
+            \gamma_1(R,\theta) = \frac{\Sigma_2}{2} - \int_R^{\infty} \mathrm{d} R'
+            \frac{\Sigma_2(R')}{R'}
+
+        where :math:`\Sigma_0,\; \Sigma_2` come from multipole coefficients of :math:`\Sigma(R,
+        \theta)`:
+
+        .. math::
+            \Sigma_0(R) = \Sigma_{\rm sph}(R)\left[1 + \dfrac{\epsilon^2}{2} \left(\eta(R) +
+            \dfrac{\eta(R)^2}{2} + \frac{1}{2}\dfrac{\mathrm{d} \eta(R)}{\mathrm{d}\ln R} \right)
+            \right] + \mathcal{O}( \epsilon^3 )
+
+        .. math::
+            \Sigma_2(R) = - \epsilon \eta(R) \Sigma_{\rm sph} (R) \cos({2\theta_0}) + \mathcal{O}(
+            \epsilon^3 ),
+
+
+        with:
+
+        .. math::
+            \eta(R) = \frac{\mathrm{d} \ln{\Sigma_{\rm sph}(R)}} {\mathrm{d} \ln{R}},\;
+            \epsilon = \dfrac{1-q}{1+q}
+
+        and :math:`\theta_0` being the polar angle of the major axis (:math:`\theta_0=0` when we
+        align the major axis to the :math:`x`-axis).
+
+        Parameters
+        ----------
+        r_proj: array
+            Projected radial position from the cluster center in :math:`M\!pc`.
+        z_cl: float
+            Redshift of lens cluster
+        ell: float
+            ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
+            q=b/a (Ratio of major axis to the minor axis lengths)
+        term: str
+            The expansion term wanted.
+                * 'mono': The monopole term with the ellipticity corrections applied. This will
+                    give the usual excess surface density but for a triaxial halo.
+                * 'quad_4theta': The 4theta component of the quadrupole term.
+                * 'quad_const': The constant component of the quadrupole term.
+        n_grid: int
+            Grid steps for gradient calculations.
+
+        Returns
+        -------
+        numpy.ndarry, float
+            Requested component of the excess surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+        """
+
+        if self.validate_input:
+            validate_argument(locals(), "r_proj", "float_array", argmin=0)
+            validate_argument(locals(), "z_cl", float, argmin=0)
+            validate_argument(locals(), "ell", float, argmin=0, argmax=1)
+            validate_argument(locals(), "term", str)
+            validate_argument(locals(), "n_grid", int, argmin=2)
+            validate_argument(locals(), "r_mis", float, argmin=0, eqmin=True, none_ok=True)
+
+        if self.halo_profile_model == "einasto" and verbose:
+            print(f"Einasto alpha = {self._get_einasto_alpha(z_cl=z_cl)}")
+
+        if r_mis is not None:
+
+            def surface_density_func(r_proj, z_cl):
+                return self._eval_surface_density_miscentered(
+                    r_proj=r_proj, z_cl=z_cl, r_mis=r_mis, mis_from_backend=mis_from_backend
+                )
+
+        else:
+            surface_density_func = self._eval_surface_density
+
+        if self.backend not in ("ccl", "nc"):
+            raise NotImplementedError(
+                f"Triaxial-4theta term not currently supported with the {self.backend} backend. "
+                "Use the CCL or NumCosmo backend instead"
+            )
+
+        return self._eval_excess_surface_density_triaxial(
+            surface_density_func, r_proj, z_cl, ell, term, n_grid
+        )
+
     def eval_tangential_shear(self, r_proj, z_cl, z_src, z_src_info="discrete", verbose=False):
         r"""Computes the tangential shear
 
diff --git a/clmm/theory/triaxiality.py b/clmm/theory/triaxiality.py
new file mode 100644
index 000000000..e239b1277
--- /dev/null
+++ b/clmm/theory/triaxiality.py
@@ -0,0 +1,250 @@
+"""@file triaxialit.py
+Model functions with triaxiality
+"""
+
+import numpy as np
+from scipy.interpolate import InterpolatedUnivariateSpline
+
+
+def _d_dlnr(r_grid, f_grid):
+    """Derivate of function with respect to ln(R)
+
+    Parameters
+    ----------
+    r_grid: np.array
+        Radius value in a grid
+    f_grid: np.array
+        Corresponding value of the function in the grid
+
+    Returns
+    -------
+    np.array
+        Derivative of function with respect to ln(R).
+    """
+    return r_grid * np.gradient(f_grid, r_grid)
+
+
+def _value_from_grid(r_grid, f_grid, r_value):
+    """Get value of tabulated function at given value
+
+    Parameters
+    ----------
+    r_grid: np.array
+        Radius value in a grid
+    f_grid: np.array
+        Corresponding value of the function in the grid
+    r_value: float
+        Value where to evaluate function
+
+    Returns
+    -------
+    float
+        Value of function at r_value.
+    """
+    return InterpolatedUnivariateSpline(r_grid, f_grid, k=5)(r_value)
+
+
+def _integrate_grid(r_grid, f_grid, r_limits):
+    """Integrate tabulated function.
+
+    Parameters
+    ----------
+    r_grid: np.array
+        Radius value in a grid
+    f_grid: np.array
+        Corresponding value of the function in the grid
+    r_limits: tuple
+        Limits of integration.
+
+    Returns
+    -------
+    float
+        Integrated value of function.
+    """
+    return np.vectorize(InterpolatedUnivariateSpline(r_grid, f_grid, k=5).integral)(*r_limits)
+
+
+def _sigma0_correction(ell, sigma_sph, eta, deta_dlnr):
+    """
+    Zero-th order correction for ellipticity.
+
+    Parameters
+    ----------
+    ell: float
+        ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
+    sigma_sph: array
+        Surface density assuming the halo is spherical.
+    eta: array
+        Derivative of sigma_sph with radius: dln(sigma_sph)/dln(R)
+    deta_dlnr: array
+        Derivative of eta with radius: dln(eta)/dln(R)
+
+    Returns
+    -------
+    float, np.array
+        Zero-th order correction for ellipticity.
+    """
+    return sigma_sph * (0.5 * ell**2 * (eta + 0.5 * eta**2 + 0.5 * deta_dlnr))
+
+
+def excess_surface_density_mono_correction(surface_density_func, r_proj, z_cl, ell, n_grid=10000):
+    r"""
+    Compute the ellipticity correction for the monopole term:
+
+    .. math::
+        \Sigma_0(R) - \Sigma_{\rm sph}(R) = \Sigma_{\rm sph}(R)\left[\dfrac{\epsilon^2}{2}
+        \left(\eta(R) + \dfrac{\eta(R)^2}{2} + \frac{1}{2}\dfrac{\mathrm{d} \eta(R)}{\mathrm{d}\ln
+        R} \right) \right] + \mathcal{O}( \epsilon^3 )
+
+
+    Parameters
+    ----------
+    surface_density_func: function
+        Function that computes the surface density, must take (r_proj, z_cl) as input.
+    r_proj: array
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    z_cl: float
+        Redshift of lens cluster
+    ell: float
+        ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
+        q=b/a (Ratio of major axis to the minor axis lengths)
+    n_grid: int
+        Grid steps for gradient calculations.
+
+    Returns
+    -------
+    array
+        The ellipticity correction for the monopole term.
+    """
+    # Compute grid for integration and interpolation
+    r_grid = np.logspace(-3, np.log10(3 * np.max(r_proj)), n_grid)
+    sigma_sph_grid = surface_density_func(r_grid, z_cl)
+    eta_grid = _d_dlnr(r_grid, np.log(sigma_sph_grid))
+    deta_dlnr_grid = _d_dlnr(r_grid, eta_grid)
+
+    # Compute equation terms
+    integral_term = _integrate_grid(
+        r_grid,
+        r_grid * _sigma0_correction(ell, sigma_sph_grid, eta_grid, deta_dlnr_grid),
+        (0, r_proj),
+    )
+    sigma_correction = _sigma0_correction(
+        ell,
+        sigma_sph=surface_density_func(r_proj, z_cl),
+        eta=_value_from_grid(r_grid, eta_grid, r_proj),
+        deta_dlnr=_value_from_grid(r_grid, deta_dlnr_grid, r_proj),
+    )
+
+    return (2 * integral_term / r_proj**2) - sigma_correction
+
+
+def excess_surface_density_quad_4theta(surface_density_func, r_proj, z_cl, ell, n_grid=10000):
+    r"""
+    Compute the 4theta component of the quadrupole term:
+
+    .. math::
+        \Delta\Sigma^{4\theta} (R) \equiv \frac{\Sigma_{\rm crit}}{2\pi} \int \mathrm{d}\theta
+        \left( \gamma_1(R,\theta)\cos{4\theta} + \gamma_2(R,\theta) \sin{4\theta} \right) \\ =
+        \frac{\Sigma_2}{2} - \frac{3}{R^4}\int_0^R \mathrm{d}R' R'^3 \Sigma_2 (R')
+
+    where :math:`\Sigma_2` come from multipole coefficients of :math:`\Sigma(R, \theta)`:
+
+    .. math::
+        \Sigma_2(R) = - \epsilon \eta(R) \Sigma_{\rm sph} (R) \cos({2\theta_0}) + \mathcal{O}(
+        \epsilon^3 ),
+
+    with:
+
+    .. math::
+        \eta(R) = \frac{\mathrm{d} \ln{\Sigma_{\rm sph}(R)}} {\mathrm{d} \ln{R}},\;
+        \epsilon = \dfrac{1-q}{1+q}
+
+    and :math:`\theta_0` being the polar angle of the major axis (:math:`\theta_0=0` when we
+    align the major axis to the :math:`x`-axis).
+
+    Parameters
+    ----------
+    surface_density_func: function
+        Function that computes the surface density, must take (r_proj, z_cl) as input.
+    r_proj: array
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    z_cl: float
+        Redshift of lens cluster
+    ell: float
+        ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
+        q=b/a (Ratio of major axis to the minor axis lengths)
+    n_grid: int
+        Grid steps for gradient calculations.
+
+    Returns
+    -------
+    array
+        The 4theta component of the quadrupole term.
+    """
+    # Compute grid for integration and interpolation
+    r_grid = np.logspace(-3, np.log10(3 * np.max(r_proj)), n_grid)
+    sigma_sph_grid = surface_density_func(r_grid, z_cl)
+    eta_grid = _d_dlnr(r_grid, np.log(sigma_sph_grid))
+
+    # Compute equation terms
+    integral_term = _integrate_grid(r_grid, r_grid**3 * sigma_sph_grid * eta_grid, (0, r_proj))
+    sigma_sph = surface_density_func(r_proj, z_cl)
+    eta = _value_from_grid(r_grid, eta_grid, r_proj)
+
+    return ell * (3 * integral_term / r_proj**4 - 0.5 * sigma_sph * eta)
+
+
+def excess_surface_density_quad_const(surface_density_func, r_proj, z_cl, ell, n_grid=10000):
+    r"""
+    Compute the constant component of the quadrupole term:
+
+    .. math::
+        \Delta\Sigma^{\rm const} (R) \equiv \frac{\Sigma_{\rm crit}}{2\pi} \int \mathrm{d}\theta
+        \gamma_1(R,\theta) = \frac{\Sigma_2}{2} - \int_R^{\infty} \mathrm{d} R'
+        \frac{\Sigma_2(R')}{R'}
+
+    where :math:`\Sigma_2` come from multipole coefficients of :math:`\Sigma(R, \theta)`:
+
+    .. math::
+        \Sigma_2(R) = - \epsilon \eta(R) \Sigma_{\rm sph} (R) \cos({2\theta_0}) + \mathcal{O}(
+        \epsilon^3 ),
+
+    with:
+
+    .. math::
+        \eta(R) = \frac{\mathrm{d} \ln{\Sigma_{\rm sph}(R)}} {\mathrm{d} \ln{R}},\;
+        \epsilon = \dfrac{1-q}{1+q}
+
+    and :math:`\theta_0` being the polar angle of the major axis (:math:`\theta_0=0` when we
+    align the major axis to the :math:`x`-axis).
+
+    Parameters
+    ----------
+    surface_density_func: function
+        Function that computes the surface density, must take (r_proj, z_cl) as input.
+    r_proj: array
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    z_cl: float
+        Redshift of lens cluster
+    ell: float
+        ellipticity of halo defined by e = (1-q)/(1+q), q is the axis ratio.
+        q=b/a (Ratio of major axis to the minor axis lengths)
+    n_grid: int
+        Grid steps for gradient calculations.
+
+    Returns
+    -------
+    array
+        The constant component of the quadrupole term.
+    """
+    # Compute grid for integration and interpolation
+    r_grid = np.logspace(-3, np.log10(3 * np.max(r_proj)), n_grid)
+    sigma_sph_grid = surface_density_func(r_grid, z_cl)
+    eta_grid = _d_dlnr(r_grid, np.log(sigma_sph_grid))
+
+    # Compute equation terms
+    integral_term = _integrate_grid(r_grid, sigma_sph_grid * eta_grid / r_grid, (r_proj, np.inf))
+    sigma_sph = surface_density_func(r_proj, z_cl)
+    eta = _value_from_grid(r_grid, eta_grid, r_proj)
+
+    return ell * (integral_term - 0.5 * sigma_sph * eta)
diff --git a/clmm/utils/validation.py b/clmm/utils/validation.py
index e30d6c4a9..0b0c75d6f 100644
--- a/clmm/utils/validation.py
+++ b/clmm/utils/validation.py
@@ -212,7 +212,7 @@ def _validate_dec(loc, dec_name, is_array):
 
 
 def _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c):
-    r""" "Validate the compatibility between is_deltasigma and sigma_c arguments.
+    r"""Validate the compatibility between is_deltasigma and sigma_c arguments.
 
 
     Parameters
diff --git a/examples/demo_dataops_functionality.ipynb b/examples/demo_dataops_functionality.ipynb
index a985af223..99b1d7d33 100644
--- a/examples/demo_dataops_functionality.ipynb
+++ b/examples/demo_dataops_functionality.ipynb
@@ -157,7 +157,11 @@
    "outputs": [],
    "source": [
     "cl = GalaxyCluster(\n",
-    "    cluster_id, cluster_ra, cluster_dec, cluster_z, noisy_data_z, coordinate_system=\"euclidean\"\n",
+    "    cluster_id,\n",
+    "    cluster_ra,\n",
+    "    cluster_dec,\n",
+    "    cluster_z,\n",
+    "    noisy_data_z,\n",
     ")"
    ]
   },
@@ -279,6 +283,52 @@
     ")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It is also possible to compute the quadrupole compenents ($\\Delta\\Sigma^{\\rm const}$, $\\Delta\\Sigma^{4\\theta}$), where the approach of\n",
+    "[Shin et al. 2018](https://doi.org/10.1093/mnras/stx3366) is used.\n",
+    "In this case `compute_tangential_and_cross_components` returns (`theta`, `e_t`, `e_x`, `e_quad_4theta`, `e_quad_const`):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "theta_and_comps_quad = compute_tangential_and_cross_components(\n",
+    "    ra_lens=cl.ra,\n",
+    "    dec_lens=cl.dec,\n",
+    "    ra_source=cl.galcat[\"ra\"],\n",
+    "    dec_source=cl.galcat[\"dec\"],\n",
+    "    shear1=cl.galcat[\"e1\"],\n",
+    "    shear2=cl.galcat[\"e2\"],\n",
+    "    coordinate_system=\"euclidean\",\n",
+    "    include_quadrupole=True,\n",
+    "    phi_major=1.0,  # orientation of major axis\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can compare that the first three components are the same"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "(theta_and_comps_quad[0] == theta).all(), (theta_and_comps_quad[1] == e_t).all(), (\n",
+    "    theta_and_comps_quad[2] == e_x\n",
+    ").all()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -323,6 +373,27 @@
     "cl.galcat[\"e_tan\", \"e_cross\"].pprint(max_width=-1)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It is also possible to also compute the quadrupole components"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cl.include_quadrupole = True\n",
+    "cl.set_phi_major(1.0)  # orientation of major axis\n",
+    "cl.compute_tangential_and_cross_components(add=True)\n",
+    "# With the option add the cl object has theta, et and ex new columns\n",
+    "# (default: takes in columns named 'e1' and 'e2' and save the results in 'et' and 'ex')\n",
+    "cl.galcat[\"et\", \"ex\", \"e_quad_4theta\", \"e_quad_const\"].pprint(max_width=-1)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -336,14 +407,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "f, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
+    "f, ax = plt.subplots(2, 2, figsize=(10, 10))\n",
+    "\n",
+    "ax[0, 0].hist(cl.galcat[\"et\"], bins=50)\n",
+    "ax[0, 0].set_xlabel(\"$\\\\epsilon_t$\", fontsize=\"xx-large\")\n",
     "\n",
-    "ax[0].hist(cl.galcat[\"et\"], bins=50)\n",
-    "ax[0].set_xlabel(\"$\\\\epsilon_t$\", fontsize=\"xx-large\")\n",
+    "ax[0, 1].hist(cl.galcat[\"ex\"], bins=50)\n",
+    "ax[0, 1].set_xlabel(\"$\\\\epsilon_x$\", fontsize=\"xx-large\")\n",
+    "ax[0, 1].set_yscale(\"log\")\n",
     "\n",
-    "ax[1].hist(cl.galcat[\"ex\"], bins=50)\n",
-    "ax[1].set_xlabel(\"$\\\\epsilon_x$\", fontsize=\"xx-large\")\n",
-    "ax[1].set_yscale(\"log\")"
+    "ax[1, 0].hist(cl.galcat[\"e_quad_4theta\"], bins=50)\n",
+    "ax[1, 0].set_xlabel(\"$\\\\epsilon_{4\\\\theta}$\", fontsize=\"xx-large\")\n",
+    "\n",
+    "ax[1, 1].hist(cl.galcat[\"e_quad_const\"], bins=50)\n",
+    "ax[1, 1].set_xlabel(\"$\\\\epsilon_{\\\\rm const}$\", fontsize=\"xx-large\")"
    ]
   },
   {
@@ -378,7 +455,7 @@
     "    z_lens=cl.z,\n",
     ")\n",
     "# profiles.pprint(max_width=-1)\n",
-    "profiles.show_in_notebook()"
+    "profiles  # .show_in_notebook()"
    ]
   },
   {
@@ -404,7 +481,7 @@
    "source": [
     "cl.make_radial_profile(\"kpc\", cosmo=cosmo)\n",
     "# cl.profile.pprint(max_width=-1)\n",
-    "cl.profile.show_in_notebook()"
+    "cl.profile  # .show_in_notebook()"
    ]
   },
   {
@@ -497,8 +574,9 @@
     "new_bins = make_bins(1, 6, nbins=20, method=\"evenlog10width\")\n",
     "\n",
     "new_profiles = cl.make_radial_profile(\"Mpc\", bins=new_bins, cosmo=cosmo)\n",
-    "fig1, ax1 = cl.plot_profiles()\n",
-    "ax1.set_xscale(\"log\")"
+    "(fig1, ax1), (fig2, ax2) = cl.plot_profiles()\n",
+    "ax1.set_xscale(\"log\")\n",
+    "ax2.set_xscale(\"log\")"
    ]
   },
   {
@@ -572,7 +650,7 @@
    "source": [
     "# Here the list of galaxy IDs that are in the first bin of the tangential shear profile\n",
     "gal_list = cl.profile[\"gal_id\"][0]\n",
-    "print(gal_list)"
+    "print(np.array(gal_list, dtype=int))"
    ]
   },
   {
@@ -598,7 +676,7 @@
     "    cross_component_out=\"g_cross\",\n",
     ")\n",
     "# cl.profile.pprint(max_width=-1)\n",
-    "cl.profile.show_in_notebook()"
+    "cl.profile  # .show_in_notebook()"
    ]
   },
   {
@@ -624,7 +702,7 @@
     "    table_name=\"reduced_shear_profile\",\n",
     ")\n",
     "# cl.reduced_shear_profile.pprint(max_width=-1)\n",
-    "cl.reduced_shear_profile.show_in_notebook()"
+    "cl.reduced_shear_profile  # .show_in_notebook()"
    ]
   },
   {
@@ -647,6 +725,8 @@
     "    shape_component2=\"e2\",\n",
     "    tan_component=\"DeltaSigma_tan\",\n",
     "    cross_component=\"DeltaSigma_cross\",\n",
+    "    quad_4theta_component=\"DeltaSigma_4theta\",\n",
+    "    quad_const_component=\"DeltaSigma_const\",\n",
     "    add=True,\n",
     "    cosmo=cosmo,\n",
     "    is_deltasigma=True,\n",
@@ -710,7 +790,7 @@
    "outputs": [],
    "source": [
     "\"\"\"\n",
-    "cl.make_radial_profile(\"Mpc\", cosmo=cosmo, \n",
+    "cl.make_radial_profile(\"Mpc\", cosmo=cosmo,\n",
     "                       tan_component_in='DeltaSigma_tan', cross_component_in='DeltaSigma_cross',\n",
     "                       tan_component_out='DeltaSigma_tan', cross_component_out='DeltaSigma_cross',\n",
     "                       table_name='DeltaSigma_profile').pprint(max_width=-1)\n",
@@ -723,6 +803,10 @@
     "    cross_component_in=\"DeltaSigma_cross\",\n",
     "    tan_component_out=\"DeltaSigma_tan\",\n",
     "    cross_component_out=\"DeltaSigma_cross\",\n",
+    "    quad_4theta_component_in=\"DeltaSigma_4theta\",\n",
+    "    quad_const_component_in=\"DeltaSigma_const\",\n",
+    "    quad_4theta_component_out=\"DeltaSigma_4theta\",\n",
+    "    quad_const_component_out=\"DeltaSigma_const\",\n",
     "    table_name=\"DeltaSigma_profile\",\n",
     "    use_weights=True,\n",
     ")"
@@ -735,7 +819,7 @@
    "outputs": [],
    "source": [
     "# cl.DeltaSigma_profile.pprint(max_width=-1)\n",
-    "cl.DeltaSigma_profile.show_in_notebook()"
+    "cl.DeltaSigma_profile  # .show_in_notebook()"
    ]
   },
   {
@@ -752,7 +836,7 @@
    "outputs": [],
    "source": [
     "avg_profile_noweights = make_radial_profile(\n",
-    "    [cl.galcat[\"DeltaSigma_tan\"], cl.galcat[\"DeltaSigma_cross\"]],\n",
+    "    [cl.galcat[\"DeltaSigma_tan\"], cl.galcat[\"DeltaSigma_cross\"], cl.galcat[\"DeltaSigma_4theta\"]],\n",
     "    cl.galcat[\"theta\"],\n",
     "    z_lens=cl.z,\n",
     "    bin_units=\"Mpc\",\n",
@@ -790,11 +874,39 @@
     "plt.legend()\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.errorbar(\n",
+    "    cl.DeltaSigma_profile[\"radius\"],\n",
+    "    cl.DeltaSigma_profile[\"DeltaSigma_4theta\"],\n",
+    "    cl.DeltaSigma_profile[\"DeltaSigma_4theta_err\"],\n",
+    "    marker=\"o\",\n",
+    "    label=\"using weights\",\n",
+    ")\n",
+    "plt.errorbar(\n",
+    "    avg_profile_noweights[\"radius\"],\n",
+    "    avg_profile_noweights[\"p_2\"],\n",
+    "    avg_profile_noweights[\"p_2_err\"],\n",
+    "    marker=\"x\",\n",
+    "    label=\"no weights\",\n",
+    ")\n",
+    "\n",
+    "plt.title(\"DeltaSigma profile\")\n",
+    "plt.xlabel(\"Radius [Mpc]\")\n",
+    "plt.ylabel(\"$\\Delta\\Sigma^{4\\\\theta} [M_\\odot\\; Mpc^{-2}]$\")\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -808,7 +920,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.9"
+   "version": "3.13.7"
   }
  },
  "nbformat": 4,
diff --git a/examples/demo_theory_functionality.ipynb b/examples/demo_theory_functionality.ipynb
index 69146545b..6493ee4df 100644
--- a/examples/demo_theory_functionality.ipynb
+++ b/examples/demo_theory_functionality.ipynb
@@ -200,6 +200,83 @@
     ")"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "DeltaSigma_triax_mono = m.compute_excess_surface_density_triaxial(\n",
+    "    r3d,\n",
+    "    cluster_mass,\n",
+    "    cluster_concentration,\n",
+    "    z_cl,\n",
+    "    ell=0.7,\n",
+    "    cosmo=cosmo,\n",
+    "    term=\"mono\",\n",
+    "    delta_mdef=mass_Delta,\n",
+    "    halo_profile_model=density_profile_parametrization,\n",
+    ")\n",
+    "DeltaSigma_triax_quad_4theta = m.compute_excess_surface_density_triaxial(\n",
+    "    r3d,\n",
+    "    cluster_mass,\n",
+    "    cluster_concentration,\n",
+    "    z_cl,\n",
+    "    ell=0.7,\n",
+    "    cosmo=cosmo,\n",
+    "    term=\"quad_4theta\",\n",
+    "    delta_mdef=mass_Delta,\n",
+    "    halo_profile_model=density_profile_parametrization,\n",
+    ")\n",
+    "DeltaSigma_triax_quad_const = m.compute_excess_surface_density_triaxial(\n",
+    "    r3d,\n",
+    "    cluster_mass,\n",
+    "    cluster_concentration,\n",
+    "    z_cl,\n",
+    "    ell=0.7,\n",
+    "    cosmo=cosmo,\n",
+    "    term=\"quad_const\",\n",
+    "    delta_mdef=mass_Delta,\n",
+    "    halo_profile_model=density_profile_parametrization,\n",
+    ")\n",
+    "DeltaSigma_triax_mono_mis = m.compute_excess_surface_density_triaxial(\n",
+    "    r3d,\n",
+    "    cluster_mass,\n",
+    "    cluster_concentration,\n",
+    "    z_cl,\n",
+    "    ell=0.7,\n",
+    "    cosmo=cosmo,\n",
+    "    term=\"mono\",\n",
+    "    delta_mdef=mass_Delta,\n",
+    "    halo_profile_model=density_profile_parametrization,\n",
+    "    r_mis=0.2,\n",
+    ")\n",
+    "DeltaSigma_triax_quad_4theta_mis = m.compute_excess_surface_density_triaxial(\n",
+    "    r3d,\n",
+    "    cluster_mass,\n",
+    "    cluster_concentration,\n",
+    "    z_cl,\n",
+    "    ell=0.7,\n",
+    "    cosmo=cosmo,\n",
+    "    term=\"quad_4theta\",\n",
+    "    delta_mdef=mass_Delta,\n",
+    "    halo_profile_model=density_profile_parametrization,\n",
+    "    r_mis=0.2,\n",
+    ")\n",
+    "DeltaSigma_triax_quad_const_mis = m.compute_excess_surface_density_triaxial(\n",
+    "    r3d,\n",
+    "    cluster_mass,\n",
+    "    cluster_concentration,\n",
+    "    z_cl,\n",
+    "    ell=0.7,\n",
+    "    cosmo=cosmo,\n",
+    "    term=\"quad_const\",\n",
+    "    delta_mdef=mass_Delta,\n",
+    "    halo_profile_model=density_profile_parametrization,\n",
+    "    r_mis=0.2,\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -349,8 +426,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def plot_profile(r, profile_vals, profile_label=\"rho\", label=\"\"):\n",
-    "    plt.loglog(r, profile_vals, label=label)\n",
+    "def plot_profile(r, profile_vals, profile_label=\"rho\", **kwargs):\n",
+    "    plt.loglog(r, profile_vals, **kwargs)\n",
     "    plt.xlabel(\"r [Mpc]\", fontsize=\"xx-large\")\n",
     "    plt.ylabel(profile_label, fontsize=\"xx-large\")"
    ]
@@ -388,6 +465,45 @@
     "plt.legend()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_profile(r3d, DeltaSigma_triax_mono, None, label=\"$\\\\Delta\\\\Sigma^{\\\\rm mono}$\")\n",
+    "plot_profile(r3d, DeltaSigma_triax_quad_4theta, None, label=\"$\\\\Delta\\\\Sigma^{4\\\\theta}$\")\n",
+    "plot_profile(r3d, DeltaSigma_triax_quad_const, None, label=\"$\\\\Delta\\\\Sigma^{\\\\rm const.}$\")\n",
+    "\n",
+    "plot_profile(\n",
+    "    r3d,\n",
+    "    DeltaSigma_triax_mono_mis,\n",
+    "    None,\n",
+    "    label=\"$\\\\Delta\\\\Sigma^{\\\\rm mono}_{\\\\rm miscent.}$\",\n",
+    "    ls=\"--\",\n",
+    "    color=\"C0\",\n",
+    ")\n",
+    "plot_profile(\n",
+    "    r3d,\n",
+    "    DeltaSigma_triax_quad_4theta_mis,\n",
+    "    None,\n",
+    "    label=\"$\\\\Delta\\\\Sigma^{4\\\\theta}_{\\\\rm miscent.}$\",\n",
+    "    ls=\"--\",\n",
+    "    color=\"C1\",\n",
+    ")\n",
+    "plot_profile(\n",
+    "    r3d,\n",
+    "    DeltaSigma_triax_quad_const_mis,\n",
+    "    \"$\\\\Delta\\\\Sigma_{\\\\rm triaxial}$\",\n",
+    "    label=\"$\\\\Delta\\\\Sigma^{\\\\rm const.}_{\\\\rm miscent.}$\",\n",
+    "    ls=\"--\",\n",
+    "    color=\"C2\",\n",
+    ")\n",
+    "plt.yscale(\"linear\")\n",
+    "plt.ylim(-0.8e15, 0.8e15)\n",
+    "plt.legend(ncol=2)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -571,9 +687,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "wrk"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -585,7 +701,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.9"
+   "version": "3.13.7"
   }
  },
  "nbformat": 4,
diff --git a/examples/demo_theory_functionality_oo.ipynb b/examples/demo_theory_functionality_oo.ipynb
index 56197d8b7..89ee924b0 100644
--- a/examples/demo_theory_functionality_oo.ipynb
+++ b/examples/demo_theory_functionality_oo.ipynb
@@ -136,6 +136,25 @@
     "gammat = moo.eval_tangential_shear(r3d, z_cl, z_src)\n",
     "kappa = moo.eval_convergence(r3d, z_cl, z_src)\n",
     "\n",
+    "# Triaxial DeltaSigam\n",
+    "DeltaSigma_triax_mono = moo.eval_excess_surface_density_triaxial(r3d, z_cl, ell=0.7, term=\"mono\")\n",
+    "DeltaSigma_triax_quad_4theta = moo.eval_excess_surface_density_triaxial(\n",
+    "    r3d, z_cl, ell=0.7, term=\"quad_4theta\"\n",
+    ")\n",
+    "DeltaSigma_triax_quad_const = moo.eval_excess_surface_density_triaxial(\n",
+    "    r3d, z_cl, ell=0.7, term=\"quad_const\"\n",
+    ")\n",
+    "DeltaSigma_triax_mono_mis = moo.eval_excess_surface_density_triaxial(\n",
+    "    r3d, z_cl, ell=0.7, term=\"mono\", r_mis=0.2\n",
+    ")\n",
+    "DeltaSigma_triax_quad_4theta_mis = moo.eval_excess_surface_density_triaxial(\n",
+    "    r3d, z_cl, ell=0.7, term=\"quad_4theta\", r_mis=0.2\n",
+    ")\n",
+    "DeltaSigma_triax_quad_const_mis = moo.eval_excess_surface_density_triaxial(\n",
+    "    r3d, z_cl, ell=0.7, term=\"quad_const\", r_mis=0.2\n",
+    ")\n",
+    "\n",
+    "\n",
     "gt = moo.eval_reduced_tangential_shear(r3d, z_cl, z_src)\n",
     "# Lensing quantities assuming sources follow a given redshift distribution.\n",
     "\n",
@@ -170,8 +189,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def plot_profile(r, profile_vals, profile_label=\"rho\", label=\"\"):\n",
-    "    plt.loglog(r, profile_vals, label=label)\n",
+    "def plot_profile(r, profile_vals, profile_label=\"rho\", **kwargs):\n",
+    "    plt.loglog(r, profile_vals, **kwargs)\n",
     "    plt.xlabel(\"r [Mpc]\", fontsize=\"xx-large\")\n",
     "    plt.ylabel(profile_label, fontsize=\"xx-large\")"
    ]
@@ -212,6 +231,45 @@
     "plt.savefig(\"miscentering.png\")"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_profile(r3d, DeltaSigma_triax_mono, None, label=\"$\\\\Delta\\\\Sigma^{\\\\rm mono}$\")\n",
+    "plot_profile(r3d, DeltaSigma_triax_quad_4theta, None, label=\"$\\\\Delta\\\\Sigma^{4\\\\theta}$\")\n",
+    "plot_profile(r3d, DeltaSigma_triax_quad_const, None, label=\"$\\\\Delta\\\\Sigma^{\\\\rm const.}$\")\n",
+    "\n",
+    "plot_profile(\n",
+    "    r3d,\n",
+    "    DeltaSigma_triax_mono_mis,\n",
+    "    None,\n",
+    "    label=\"$\\\\Delta\\\\Sigma^{\\\\rm mono}_{\\\\rm miscent.}$\",\n",
+    "    ls=\"--\",\n",
+    "    color=\"C0\",\n",
+    ")\n",
+    "plot_profile(\n",
+    "    r3d,\n",
+    "    DeltaSigma_triax_quad_4theta_mis,\n",
+    "    None,\n",
+    "    label=\"$\\\\Delta\\\\Sigma^{4\\\\theta}_{\\\\rm miscent.}$\",\n",
+    "    ls=\"--\",\n",
+    "    color=\"C1\",\n",
+    ")\n",
+    "plot_profile(\n",
+    "    r3d,\n",
+    "    DeltaSigma_triax_quad_const_mis,\n",
+    "    \"$\\\\Delta\\\\Sigma_{\\\\rm triaxial}$\",\n",
+    "    label=\"$\\\\Delta\\\\Sigma^{\\\\rm const.}_{\\\\rm miscent.}$\",\n",
+    "    ls=\"--\",\n",
+    "    color=\"C2\",\n",
+    ")\n",
+    "plt.yscale(\"linear\")\n",
+    "plt.ylim(-0.8e15, 0.8e15)\n",
+    "plt.legend(ncol=2)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -369,9 +427,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "wrk",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "wrk"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -383,7 +441,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.13.7"
   }
  },
  "nbformat": 4,
diff --git a/tests/test_clusterensemble.py b/tests/test_clusterensemble.py
index 63c3d3958..ed1397371 100644
--- a/tests/test_clusterensemble.py
+++ b/tests/test_clusterensemble.py
@@ -1,6 +1,7 @@
 """
 tests for clusterensemble.py
 """
+
 import os
 from numpy.testing import assert_raises, assert_equal, assert_allclose, assert_array_equal
 import clmm
@@ -39,9 +40,17 @@ def test_cluster_ensemble():
     cluster = clmm.GalaxyCluster(
         unique_id="test", ra=ra_lens, dec=dec_lens, z=z_lens, galcat=galcat
     )
+    cluster_quad = clmm.GalaxyCluster(
+        unique_id="test", ra=ra_lens, dec=dec_lens, z=z_lens, galcat=galcat, include_quadrupole=True
+    )
     cluster.compute_tangential_and_cross_components()
+    cluster_quad.set_phi_major(
+        info_mem=(np.array([119.9, 120.1]), np.array([41.9, 42.1]), np.array([1.0, 1.0]))
+    )
+    cluster_quad.compute_tangential_and_cross_components()
     bins = bins_radians
     gc_list = [cluster]
+    gc_list_quad = [cluster_quad]
     # check empty cluster list
     ce_empty = ClusterEnsemble(
         "cluster_ensemble",
@@ -52,7 +61,20 @@ def test_cluster_ensemble():
         bin_units="radians",
         cosmo=cosmo,
     )
+    ce_empty_quad = ClusterEnsemble(
+        "cluster_ensemble",
+        tan_component_in="et",
+        cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
+        weights_in="w_ls",
+        bins=bins,
+        bin_units="radians",
+        cosmo=cosmo,
+        include_quadrupole=True,
+    )
     assert_raises(ValueError, ce_empty.make_stacked_radial_profile)
+    assert_raises(ValueError, ce_empty_quad.make_stacked_radial_profile)
     ce_empty.make_individual_radial_profile(
         cluster,
         tan_component_in="et",
@@ -62,9 +84,21 @@ def test_cluster_ensemble():
         bin_units="radians",
         cosmo=cosmo,
     )
+    ce_empty_quad.make_individual_radial_profile(
+        cluster_quad,
+        tan_component_in="et",
+        cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
+        weights_in="w_ls",
+        bins=bins,
+        bin_units="radians",
+        cosmo=cosmo,
+    )
 
     # check bad id
     assert_raises(TypeError, ClusterEnsemble, 1.3, gc_list)
+    assert_raises(TypeError, ClusterEnsemble, 1.3, gc_list_quad)
 
     # test without kwargs, args
     ce = ClusterEnsemble(
@@ -77,9 +111,23 @@ def test_cluster_ensemble():
         bin_units="radians",
         cosmo=cosmo,
     )
+    ce_quad = ClusterEnsemble(
+        "cluster_ensemble",
+        gc_list_quad,
+        tan_component_in="et",
+        cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
+        weights_in="w_ls",
+        bins=bins,
+        bin_units="radians",
+        cosmo=cosmo,
+        include_quadrupole=True,
+    )
 
     # test the lenght of the clusterensemble data attribute
     assert_equal(ce.__len__(), 1)
+    assert_equal(ce_quad.__len__(), 1)
 
     # test the lenght of the clusterensemble data attribute (after doubling the number of individual cluster)
     ce._add_values(
@@ -91,10 +139,24 @@ def test_cluster_ensemble():
         bin_units="radians",
         cosmo=cosmo,
     )
+    ce_quad._add_values(
+        [cluster_quad],
+        tan_component_in="et",
+        cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
+        weights_in="w_ls",
+        bins=bins,
+        bin_units="radians",
+        cosmo=cosmo,
+    )
     assert_equal(ce.__len__(), 2)
+    assert_equal(ce_quad.__len__(), 2)
     # test if the len of averaged profile has the lenght of binning axis
     assert_equal(len(ce.data["W_l"][0]), len(bins_radians) - 1)
     assert_equal(ce.__getitem__("gt"), ce.data["gt"])
+    assert_equal(len(ce_quad.data["W_l"][0]), len(bins_radians) - 1)
+    assert_equal(ce_quad.__getitem__("g_quad_const"), ce_quad.data["g_quad_const"])
 
 
 def test_covariance():
@@ -107,6 +169,7 @@ def test_covariance():
     sindec = [-0.1, 0, 0.1]
     cluster_dec = np.arcsin(sindec) * 180 / np.pi  # from -90 to 90 deg
     gclist = []
+    gclist_quad = []
     Rmin, Rmax = 0.3, 5  # Mpc
     thetamin, thetamax = Rmin / dz, Rmax / dz  # radians
     phi = np.pi
@@ -120,16 +183,33 @@ def test_covariance():
         bin_units="radians",
         cosmo=cosmo,
     )
+    ce_empty_quad = ClusterEnsemble(
+        "2",
+        tan_component_in="et",
+        cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
+        weights_in="w_ls",
+        bins=bins,
+        bin_units="radians",
+        cosmo=cosmo,
+        include_quadrupole=True,
+    )
 
     # check empty cluster list
     assert_raises(ValueError, ce_empty.compute_sample_covariance)
     assert_raises(ValueError, ce_empty.compute_bootstrap_covariance)
     assert_raises(ValueError, ce_empty.compute_jackknife_covariance)
+    assert_raises(ValueError, ce_empty_quad.compute_sample_covariance)
+    assert_raises(ValueError, ce_empty_quad.compute_bootstrap_covariance)
+    assert_raises(ValueError, ce_empty_quad.compute_jackknife_covariance)
 
     for i in range(n_catalogs):
         # generate random catalog
         e1, e2 = np.random.randn(ngals) * 0.001, np.random.randn(ngals) * 0.001
         et, ex = da._compute_tangential_shear(e1, e2, phi), da._compute_cross_shear(e1, e2, phi)
+        eft = da._compute_4theta_shear(e1, e2, phi)
+        ecn = e1
         z_gal = np.random.random(ngals) * (3 - 1.1) + 1.1
         id_gal = np.arange(ngals)
         theta_gal = np.linspace(0, 1, ngals) * (thetamax - thetamin) + thetamin
@@ -142,10 +222,21 @@ def test_covariance():
             "e2": e2,
             "et": et,
             "ex": ex,
+            "e_quad_4theta": eft,
+            "e_quad_const": ecn,
             "w_ls": w_ls,
         }
         cl = clmm.GalaxyCluster("mock_cluster", cluster_ra[i], cluster_dec[i], 1.0, GCData(data))
+        cl_quad = clmm.GalaxyCluster(
+            "mock cluster",
+            cluster_ra[i],
+            cluster_dec[i],
+            1.0,
+            GCData(data),
+            include_quadrupole=True,
+        )
         gclist.append(cl)
+        gclist_quad.append(cl_quad)
         ce_empty.make_individual_radial_profile(
             galaxycluster=cl,
             tan_component_in="et",
@@ -155,13 +246,26 @@ def test_covariance():
             bin_units="Mpc",
             cosmo=cosmo,
         )
+        ce_empty_quad.make_individual_radial_profile(
+            galaxycluster=cl_quad,
+            tan_component_in="et",
+            cross_component_in="ex",
+            quad_4theta_component_in="e_quad_4theta",
+            quad_const_component_in="e_quad_const",
+            weights_in="w_ls",
+            bins=bins,
+            bin_units="Mpc",
+            cosmo=cosmo,
+        )
 
     ensemble_id = 1
-    names = ["id", "ra", "dec", "z", "radius", "gt", "gx", "W_l"]
+    names = ["id", "ra", "dec", "z", "radius", "gt", "gx", "g_quad_4theta", "g_quad_const", "W_l"]
 
     # test without args, kwargs
     ce = ClusterEnsemble(ensemble_id)
+    ce_quad = ClusterEnsemble(ensemble_id, include_quadrupole=True)
     assert_raises(ValueError, ce.make_stacked_radial_profile)
+    assert_raises(ValueError, ce_quad.make_stacked_radial_profile)
 
     # test with args, kwargs
     ce = ClusterEnsemble(
@@ -174,50 +278,101 @@ def test_covariance():
         bin_units="Mpc",
         cosmo=cosmo,
     )
+    ce_quad = ClusterEnsemble(
+        ensemble_id,
+        gclist_quad,
+        tan_component_in="et",
+        cross_component_in="ex",
+        quad_4theta_component_in="e_quad_4theta",
+        quad_const_component_in="e_quad_const",
+        weights_in="w_ls",
+        bins=bins,
+        bin_units="Mpc",
+        cosmo=cosmo,
+        include_quadrupole=True,
+    )
 
     ce.make_stacked_radial_profile()
+    ce_quad.make_stacked_radial_profile()
 
     assert_raises(ValueError, ce.make_individual_radial_profile, gclist[0], bin_units="radians")
+    assert_raises(
+        ValueError, ce_quad.make_individual_radial_profile, gclist_quad[0], bin_units="radians"
+    )
+
     # test if te list object matches the calculation from the object with manually added clusters
     ce_empty.make_stacked_radial_profile()
+    ce_empty_quad.make_stacked_radial_profile()
     assert_array_equal(ce_empty.stacked_data, ce.stacked_data)
+    assert_array_equal(ce_empty_quad.stacked_data, ce_quad.stacked_data)
 
     # comparing brut force calculation for cross and tangential component
     gt_individual, gx_individual = ce.data["gt"], ce.data["gx"]
+    gft_individual, gcn_individual = ce_quad.data["g_quad_4theta"], ce_quad.data["g_quad_const"]
     Wl_individual = ce.data["W_l"]
     gt_stack = np.average(gt_individual, weights=Wl_individual, axis=0)
     gx_stack = np.average(gx_individual, weights=Wl_individual, axis=0)
+    gft_stack = np.average(gft_individual, weights=Wl_individual, axis=0)
+    gcn_stack = np.average(gcn_individual, weights=Wl_individual, axis=0)
     gt_stack_method = ce.stacked_data["gt"]
     gx_stack_method = ce.stacked_data["gx"]
+    gft_stack_method = ce_quad.stacked_data["g_quad_4theta"]
+    gcn_stack_method = ce_quad.stacked_data["g_quad_const"]
     assert_equal(gt_stack, gt_stack_method)
     assert_equal(gx_stack, gx_stack_method)
+    assert_equal(gft_stack, gft_stack_method)
+    assert_equal(gcn_stack, gcn_stack_method)
     ce.compute_sample_covariance()
+    ce_quad.compute_sample_covariance()
     ce.compute_bootstrap_covariance(n_bootstrap=3)
+    ce_quad.compute_bootstrap_covariance(n_bootstrap=3)
     ce.compute_jackknife_covariance(n_side=2)
+    ce_quad.compute_jackknife_covariance(n_side=2)
 
     # cross vs tangential covariances within a method -> shapes
     assert_equal(ce.cov["tan_sc"].shape, ce.cov["cross_sc"].shape)
     assert_equal(ce.cov["tan_bs"].shape, ce.cov["cross_bs"].shape)
     assert_equal(ce.cov["tan_jk"].shape, ce.cov["cross_jk"].shape)
+    # 4theta vs constant covariances within a method -> shapes
+    assert_equal(ce_quad.cov["quad_4theta_sc"].shape, ce_quad.cov["quad_const_sc"].shape)
+    assert_equal(ce_quad.cov["quad_4theta_bs"].shape, ce_quad.cov["quad_const_bs"].shape)
+    assert_equal(ce_quad.cov["quad_4theta_jk"].shape, ce_quad.cov["quad_const_jk"].shape)
 
     # comparing covariance estimation methods -> shapes
     assert_equal(ce.cov["tan_sc"].shape, ce.cov["tan_bs"].shape)
     assert_equal(ce.cov["tan_bs"].shape, ce.cov["tan_jk"].shape)
     assert_equal(ce.cov["tan_jk"].shape, ce.cov["tan_sc"].shape)
+    # comparing covariance estimation methods -> shapes
+    assert_equal(ce_quad.cov["quad_4theta_sc"].shape, ce_quad.cov["quad_4theta_bs"].shape)
+    assert_equal(ce_quad.cov["quad_4theta_bs"].shape, ce_quad.cov["quad_4theta_jk"].shape)
+    assert_equal(ce_quad.cov["quad_4theta_jk"].shape, ce_quad.cov["quad_4theta_sc"].shape)
 
-    # comparing brut force calculation for sample variance
+    # comparing brute force calculation for sample variance
     std_gt_stack = np.std(gt_individual, axis=0) / np.sqrt(n_catalogs - 1)
     assert_allclose(ce.cov["tan_sc"].diagonal() ** 0.5, std_gt_stack, 1e-6)
+    std_gft_stack = np.std(gft_individual, axis=0) / np.sqrt(n_catalogs - 1)
+    assert_allclose(ce_quad.cov["quad_4theta_sc"].diagonal() ** 0.5, std_gft_stack, 1e-6)
 
     # test save/load
     ce.save("ce.test.pkl")
+    ce_quad.save("ce_quad.test.pkl")
 
     ce2 = ClusterEnsemble.load("ce.test.pkl")
+    ce_quad2 = ClusterEnsemble.load("ce_quad.test.pkl")
     os.system("rm ce.test.pkl")
+    os.system("rm ce_quad.test.pkl")
+
     assert_array_equal(ce.stacked_data, ce2.stacked_data)
+    assert_array_equal(ce_quad.stacked_data, ce_quad2.stacked_data)
     assert_equal(ce.cov["tan_sc"].shape, ce2.cov["tan_sc"].shape)
     assert_equal(ce.cov["tan_bs"].shape, ce2.cov["tan_bs"].shape)
     assert_equal(ce.cov["tan_jk"].shape, ce2.cov["tan_jk"].shape)
     assert_equal(ce.cov["cross_sc"].shape, ce2.cov["cross_sc"].shape)
     assert_equal(ce.cov["cross_bs"].shape, ce2.cov["cross_bs"].shape)
     assert_equal(ce.cov["cross_jk"].shape, ce2.cov["cross_jk"].shape)
+    assert_equal(ce_quad.cov["quad_4theta_sc"].shape, ce_quad2.cov["quad_4theta_sc"].shape)
+    assert_equal(ce_quad.cov["quad_4theta_bs"].shape, ce_quad2.cov["quad_4theta_bs"].shape)
+    assert_equal(ce_quad.cov["quad_4theta_jk"].shape, ce_quad2.cov["quad_4theta_jk"].shape)
+    assert_equal(ce_quad.cov["quad_const_sc"].shape, ce_quad2.cov["quad_const_sc"].shape)
+    assert_equal(ce_quad.cov["quad_const_bs"].shape, ce_quad2.cov["quad_const_bs"].shape)
+    assert_equal(ce_quad.cov["quad_const_jk"].shape, ce_quad2.cov["quad_const_jk"].shape)
diff --git a/tests/test_dataops.py b/tests/test_dataops.py
index 9031005b8..c3741ab6e 100644
--- a/tests/test_dataops.py
+++ b/tests/test_dataops.py
@@ -55,6 +55,84 @@ def test_compute_tangential_shear():
     assert_allclose(da._compute_tangential_shear(0.0, 0.0, 0.3), 0.0, **TOLERANCE)
 
 
+def test_compute_4theta_shear():
+    """test compute quadrupole 4theta shear"""
+    shear1, shear2, phi = 0.15, 0.08, 0.52
+    expected_4theta_shear = -0.003271676989552594
+    four_theta_shear = da._compute_4theta_shear(shear1, shear2, phi)
+    assert_allclose(four_theta_shear, expected_4theta_shear)
+
+    shear1 = np.array([0.15, 0.40])
+    shear2 = np.array([0.08, 0.30])
+    phi = np.array([0.52, 1.23])
+    expected_4theta_shear = [-0.003271676989552594, -0.2111087147687582]
+    four_theta_shear = da._compute_4theta_shear(shear1, shear2, phi)
+    assert_allclose(four_theta_shear, expected_4theta_shear)
+
+    # test for reasonable values
+    assert_allclose(da._compute_4theta_shear(100.0, 0.0, 0.0), 100.0, **TOLERANCE)
+    assert_allclose(da._compute_4theta_shear(0.0, 100.0, np.pi / 8.0), 100.0, **TOLERANCE)
+    assert_allclose(da._compute_4theta_shear(0.0, 0.0, 0.3), 0.0, **TOLERANCE)
+
+
+def test_calculate_major_axis():
+    """test calculate major axis"""
+    ra_lens, dec_lens = 180.0, 0.0
+
+    ra_mem, dec_mem, weight_mem = [180.0, 180.0], [-0.5, 0.5], [1.0, 1.0]
+    expected_major_axis = np.pi / 2.0
+    assert_allclose(
+        da.calculate_major_axis(ra_lens, dec_lens, ra_mem, dec_mem, weight_mem),
+        expected_major_axis,
+        **TOLERANCE,
+    )
+
+    ra_mem, dec_mem, weight_mem = [179.5, 180.5], [0.0, 0.0], [1.0, 1.0]
+    expected_major_axis = 0.0
+    assert_allclose(
+        da.calculate_major_axis(ra_lens, dec_lens, ra_mem, dec_mem, weight_mem),
+        expected_major_axis,
+        **TOLERANCE,
+    )
+
+    ra_mem, dec_mem, weight_mem = [179.99, 180.01], [-0.01, 0.01], [1.0, 1.0]
+    expected_major_axis = -np.pi / 4.0
+    assert_allclose(
+        da.calculate_major_axis(ra_lens, dec_lens, ra_mem, dec_mem, weight_mem),
+        expected_major_axis,
+        **TOLERANCE,
+    )
+
+    ra_mem, dec_mem, weight_mem = [179.99, 180.01], [0.01, -0.01], [1.0, 1.0]
+    expected_major_axis = np.pi / 4.0
+    assert_allclose(
+        da.calculate_major_axis(ra_lens, dec_lens, ra_mem, dec_mem, weight_mem),
+        expected_major_axis,
+        **TOLERANCE,
+    )
+
+
+def test_rotate_shear():
+    """test rotate shear components"""
+    shear1, shear2, phi = 0.15, 0.08, 0.52
+    phi_major_45 = np.pi / 4.0
+    phi_major_90 = np.pi / 2.0
+    phi_major_180 = np.pi
+    expected_shear1_45, expected_shear2_45 = 0.08, -0.15
+    expected_shear1_90, expected_shear2_90 = -0.15, -0.08
+    expected_shear1_180, expected_shear2_180 = 0.15, 0.08
+
+    shear1_45, shear2_45 = da._rotate_shear(shear1, shear2, phi_major_45)
+    shear1_90, shear2_90 = da._rotate_shear(shear1, shear2, phi_major_90)
+    shear1_180, shear2_180 = da._rotate_shear(shear1, shear2, phi_major_180)
+
+    assert_allclose([shear1_45, shear2_45], [expected_shear1_45, expected_shear2_45], **TOLERANCE)
+    assert_allclose([shear1_90, shear2_90], [expected_shear1_90, expected_shear2_90], **TOLERANCE)
+    assert_allclose(
+        [shear1_180, shear2_180], [expected_shear1_180, expected_shear2_180], **TOLERANCE
+    )
+
+
 def test_compute_lensing_angles_flatsky():
     """test compute lensing angles flatsky"""
     ra_l, dec_l = 161.0, 65.0
@@ -200,6 +278,7 @@ def test_compute_tangential_and_cross_components(modeling_data):
     # Input values
     reltol = modeling_data["dataops_reltol"]
     ra_lens, dec_lens, z_lens = 120.0, 42.0, 0.5
+    phi_major = 0.0
     gals = GCData(
         {
             "ra": np.array([120.1, 119.9]),
@@ -215,20 +294,30 @@ def test_compute_tangential_and_cross_components(modeling_data):
         "angsep": np.array([0.0021745039090962414, 0.0037238407383072053]),
         "cross_shear": np.array([0.2780316984090899, 0.6398792901134982]),
         "tangential_shear": np.array([-0.22956126563459447, -0.02354769805831558]),
+        "four_theta_shear": np.array([-0.3324295, -0.43566851]),
+        "const_shear": np.array([0.2, 0.4]),
         # DeltaSigma expected values for clmm.Cosmology(H0=70.0, Omega_dm0=0.275, Omega_b0=0.025)
         "cross_DS": np.array([8.58093068e14, 1.33131522e15]),
         # [1224.3326297393244, 1899.6061989365176])*0.7*1.0e12*1.0002565513832675
         "tangential_DS": np.array([-7.08498103e14, -4.89926917e13]),
         # [-1010.889584349285, -69.9059242788237])*0.7*1.0e12*1.0002565513832675
+        "four_theta_DS": np.array([-1.02598173e15, -9.06439876e14]),
+        "const_DS": np.array([6.17262745e14, 8.32228959e14]),
     }
+
     expected_curve = {
         "angsep": np.array([0.002175111279323424171, 0.003723129781247932167]),
         "cross_shear": np.array([0.277590689496438781, 0.639929479722048944]),
         "tangential_shear": np.array([-0.23009434826803484841, -0.02214183783401518779]),
+        "four_theta_shear": np.array([-0.33189127, -0.43360245]),
+        "const_shear": np.array([0.2, 0.4]),
         # DeltaSigma expected values for clmm.Cosmology(H0=70.0, Omega_dm0=0.275, Omega_b0=0.025)
         "cross_DS": np.array([8.56731976e14, 1.33141964e15]),
         "tangential_DS": np.array([-7.10143363e14, -4.60676976e13]),
+        "four_theta_DS": np.array([-1.02432059e15, -9.02141297e14]),
+        "const_DS": np.array([6.17262745e14, 8.32228959e14]),
     }
+
     # Geometries to test
     geo_tests = [("flat", expected_flat), ("curve", expected_curve)]
     # Test domains on inputs
@@ -363,20 +452,102 @@ def test_compute_tangential_and_cross_components(modeling_data):
             angsep,
             expected["angsep"],
             **TOLERANCE,
-            err_msg="Angular Separation not correct when passing lists",
+            err_msg="Angular Separation not correct when passing arrays",
         )
         assert_allclose(
             tshear,
             expected["tangential_shear"],
             **TOLERANCE,
-            err_msg="Tangential Shear not correct when passing lists",
+            err_msg="Tangential Shear not correct when passing arrays",
         )
         assert_allclose(
             xshear,
             expected["cross_shear"],
             **TOLERANCE,
-            err_msg="Cross Shear not correct when passing lists",
+            err_msg="Cross Shear not correct when passing arrays",
+        )
+
+        ## Turn on quadrupole option
+        angsep, _, _, ftshear, cnshear = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=gals["ra"],
+            dec_source=gals["dec"],
+            shear1=gals["e1"],
+            shear2=gals["e2"],
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=0.0,
+        )
+        assert_allclose(
+            ftshear,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when passing arrays and phi_major",
+        )
+        assert_allclose(
+            cnshear,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when passing arrays and phi_major",
+        )
+        ## Turn on quadrupole option with Celestial coordinates
+        angsep, _, _, ftshear, cnshear = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=gals["ra"],
+            dec_source=gals["dec"],
+            shear1=gals["e1"],
+            shear2=-gals["e2"],
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=0.0,
+            coordinate_system="celestial",
+        )
+        assert_allclose(
+            ftshear,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when passing phi_major for celestial coord.",
+        )
+        assert_allclose(
+            cnshear,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when passing phi_major for celestial coord.",
+        )
+
+        ## Pass members info instead of major axis directly
+        angsep, _, _, ftshear, cnshear = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=gals["ra"],
+            dec_source=gals["dec"],
+            shear1=gals["e1"],
+            shear2=gals["e2"],
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=da.calculate_major_axis(
+                ra_lens,
+                dec_lens,
+                np.array([119.99, 120.01]),
+                np.array([42.0, 42.0]),
+                np.array([1.0, 1.0]),
+            ),
+        )
+        assert_allclose(
+            ftshear,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when passing arrays and mem_info",
         )
+        assert_allclose(
+            cnshear,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when passing arrays and mem_info",
+        )
+
         # Pass LISTS into function
         angsep, tshear, xshear = da.compute_tangential_and_cross_components(
             ra_lens=ra_lens,
@@ -405,6 +576,63 @@ def test_compute_tangential_and_cross_components(modeling_data):
             **TOLERANCE,
             err_msg="Cross Shear not correct when passing lists",
         )
+
+        ## Turn on quadrupole option
+        angsep, _, _, ftshear, cnshear = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=list(gals["ra"]),
+            dec_source=list(gals["dec"]),
+            shear1=list(gals["e1"]),
+            shear2=list(gals["e2"]),
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=0.0,
+        )
+        assert_allclose(
+            ftshear,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when passing lists and phi_major",
+        )
+        assert_allclose(
+            cnshear,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when passing listss and phi_major",
+        )
+
+        ## Pass members info instead of major axis directly
+        angsep, _, _, ftshear, cnshear = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=list(gals["ra"]),
+            dec_source=list(gals["dec"]),
+            shear1=list(gals["e1"]),
+            shear2=list(gals["e2"]),
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=da.calculate_major_axis(
+                ra_lens,
+                dec_lens,
+                np.array([119.99, 120.01]),
+                np.array([42.0, 42.0]),
+                np.array([1.0, 1.0]),
+            ),
+        )
+        assert_allclose(
+            ftshear,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when passing lists and mem_info",
+        )
+        assert_allclose(
+            cnshear,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when passing lists and mem_info",
+        )
+
         # Test without validation
         angsep, tshear, xshear = da.compute_tangential_and_cross_components(
             ra_lens=ra_lens,
@@ -434,6 +662,80 @@ def test_compute_tangential_and_cross_components(modeling_data):
             **TOLERANCE,
             err_msg="Cross Shear not correct when passing lists",
         )
+
+        ## Turn on quadrupole option
+        angsep, _, _, ftshear, cnshear = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=list(gals["ra"]),
+            dec_source=list(gals["dec"]),
+            shear1=list(gals["e1"]),
+            shear2=list(gals["e2"]),
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=0.0,
+            validate_input=False,
+        )
+        assert_allclose(
+            ftshear,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when passing lists and phi_major",
+        )
+        assert_allclose(
+            cnshear,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when passing listss and phi_major",
+        )
+
+        ## Pass members info instead of major axis directly
+        angsep, _, _, ftshear, cnshear = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=list(gals["ra"]),
+            dec_source=list(gals["dec"]),
+            shear1=list(gals["e1"]),
+            shear2=list(gals["e2"]),
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=da.calculate_major_axis(
+                ra_lens,
+                dec_lens,
+                np.array([119.99, 120.01]),
+                np.array([42.0, 42.0]),
+                np.array([1.0, 1.0]),
+            ),
+            validate_input=False,
+        )
+        assert_allclose(
+            ftshear,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when passing lists and mem_info",
+        )
+        assert_allclose(
+            cnshear,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when passing lists and mem_info",
+        )
+
+        # Test for ValueError if neither phi_major or mem positions are given with quadrupole
+        assert_raises(
+            ValueError,
+            da.compute_tangential_and_cross_components,
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=list(gals["ra"]),
+            dec_source=list(gals["dec"]),
+            shear1=list(gals["e1"]),
+            shear2=list(gals["e2"]),
+            geometry=geometry,
+            include_quadrupole=True,
+            validate_input=False,
+        )
+
         # Test without validation and float arguments
         angsep, tshear, xshear = da.compute_tangential_and_cross_components(
             ra_lens=ra_lens,
@@ -491,6 +793,14 @@ def test_compute_tangential_and_cross_components(modeling_data):
     cluster = clmm.GalaxyCluster(
         unique_id="blah", ra=ra_lens, dec=dec_lens, z=z_lens, galcat=gals["ra", "dec", "e1", "e2"]
     )
+    cluster_quad = clmm.GalaxyCluster(
+        unique_id="blah",
+        ra=ra_lens,
+        dec=dec_lens,
+        z=z_lens,
+        galcat=gals["ra", "dec", "e1", "e2"],
+        include_quadrupole=True,
+    )
     # Test error with bad name/missing column
     assert_raises(
         TypeError, cluster.compute_tangential_and_cross_components, shape_component1="crazy name"
@@ -518,6 +828,44 @@ def test_compute_tangential_and_cross_components(modeling_data):
             **TOLERANCE,
             err_msg="Cross Shear not correct when using cluster method",
         )
+        # include_quadrupole=True, with phi_major input
+        cluster_quad.set_phi_major(
+            info_mem=(np.array([119.99, 120.01]), np.array([42.0, 42.0]), np.array([1.0, 1.0]))
+        )
+        _, _, _, ftshear2, cnshear2 = cluster_quad.compute_tangential_and_cross_components(
+            geometry=geometry,
+        )
+        assert_allclose(
+            ftshear2,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when using cluster method w/ phi_major",
+        )
+        assert_allclose(
+            cnshear2,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when using cluster method w/ phi_major",
+        )
+        # include_quadrupole=True, with info_mem instead of phi_major input
+        cluster_quad.set_phi_major(
+            info_mem=(np.array([119.99, 120.01]), np.array([42.0, 42.0]), np.array([1.0, 1.0]))
+        )
+        _, _, _, ftshear3, cnshear3 = cluster_quad.compute_tangential_and_cross_components(
+            geometry=geometry,
+        )
+        assert_allclose(
+            ftshear3,
+            expected["four_theta_shear"],
+            **TOLERANCE,
+            err_msg="4theta Shear not correct when using cluster method w/ info_mem",
+        )
+        assert_allclose(
+            cnshear3,
+            expected["const_shear"],
+            **TOLERANCE,
+            err_msg="Constant Shear not correct when using cluster method w/ info_mem",
+        )
     # Check behaviour for the deltasigma option.
     cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.275, Omega_b0=0.025)
 
@@ -548,6 +896,18 @@ def test_compute_tangential_and_cross_components(modeling_data):
         is_deltasigma=True,
         sigma_c=None,
     )
+    # check validation between include_quadrupole and {phi_major|info_mem}
+    assert_raises(
+        ValueError,
+        da.compute_tangential_and_cross_components,
+        ra_lens=ra_lens,
+        dec_lens=dec_lens,
+        ra_source=gals["ra"],
+        dec_source=gals["dec"],
+        shear1=gals["e1"],
+        shear2=gals["e2"],
+        include_quadrupole=True,
+    )
     # test values
     for geometry, expected in geo_tests:
         angsep_DS, tDS, xDS = da.compute_tangential_and_cross_components(
@@ -568,6 +928,52 @@ def test_compute_tangential_and_cross_components(modeling_data):
             tDS, expected["tangential_DS"], reltol, err_msg="Tangential Shear not correct"
         )
         assert_allclose(xDS, expected["cross_DS"], reltol, err_msg="Cross Shear not correct")
+        ## Turn on include_quadrupole w/ phi_major
+        _, _, _, ftDS, cnDS = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=gals["ra"],
+            dec_source=gals["dec"],
+            shear1=gals["e1"],
+            shear2=gals["e2"],
+            is_deltasigma=True,
+            sigma_c=sigma_c,
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=0.0,
+        )
+        assert_allclose(
+            ftDS, expected["four_theta_DS"], 10 * reltol, err_msg="4theta shear not correct"
+        )
+        assert_allclose(
+            cnDS, expected["const_DS"], 10 * reltol, err_msg="constant shear not correct"
+        )
+        ## Turn on include_quadrupole w/ info_mem
+        _, _, _, ftDS, cnDS = da.compute_tangential_and_cross_components(
+            ra_lens=ra_lens,
+            dec_lens=dec_lens,
+            ra_source=gals["ra"],
+            dec_source=gals["dec"],
+            shear1=gals["e1"],
+            shear2=gals["e2"],
+            is_deltasigma=True,
+            sigma_c=sigma_c,
+            geometry=geometry,
+            include_quadrupole=True,
+            phi_major=da.calculate_major_axis(
+                ra_lens,
+                dec_lens,
+                np.array([119.99, 120.01]),
+                np.array([42.0, 42.0]),
+                np.array([1.0, 1.0]),
+            ),
+        )
+        assert_allclose(
+            ftDS, expected["four_theta_DS"], 10 * reltol, err_msg="4theta shear not correct"
+        )
+        assert_allclose(
+            cnDS, expected["const_DS"], 10 * reltol, err_msg="constant shear not correct"
+        )
     # Tests with the cluster object
     # cluster object missing source redshift, and function call missing cosmology
     cluster = clmm.GalaxyCluster(
@@ -582,6 +988,14 @@ def test_compute_tangential_and_cross_components(modeling_data):
         z=z_lens,
         galcat=gals["ra", "dec", "e1", "e2", "z"],
     )
+    cluster_quad = clmm.GalaxyCluster(
+        unique_id="blah",
+        ra=ra_lens,
+        dec=dec_lens,
+        z=z_lens,
+        galcat=gals["ra", "dec", "e1", "e2", "z"],
+        include_quadrupole=True,
+    )
     assert_raises(TypeError, cluster.compute_tangential_and_cross_components, is_deltasigma=True)
     # check values for DeltaSigma
     for geometry, expected in geo_tests:
@@ -606,6 +1020,48 @@ def test_compute_tangential_and_cross_components(modeling_data):
             reltol,
             err_msg="Cross Shear not correct when using cluster method",
         )
+        # Turn on include_quadrupole w/ phi_major
+        cluster_quad.set_phi_major(
+            info_mem=(np.array([119.99, 120.01]), np.array([42.0, 42.0]), np.array([1.0, 1.0]))
+        )
+        _, _, _, ftDS, cnDS = cluster_quad.compute_tangential_and_cross_components(
+            cosmo=cosmo,
+            is_deltasigma=True,
+            geometry=geometry,
+        )
+        assert_allclose(
+            ftDS,
+            expected["four_theta_DS"],
+            10 * reltol,
+            err_msg="4theta Shear not correct when using cluster method w/ phi_major",
+        )
+        assert_allclose(
+            cnDS,
+            expected["const_DS"],
+            10 * reltol,
+            err_msg="Constant Shear not correct when using cluster method w/ phi_major",
+        )
+        # Turn on include_quadrupole w/ info_mem
+        cluster_quad.set_phi_major(
+            info_mem=(np.array([119.99, 120.01]), np.array([42.0, 42.0]), np.array([1.0, 1.0]))
+        )
+        _, _, _, ftDS, cnDS = cluster_quad.compute_tangential_and_cross_components(
+            cosmo=cosmo,
+            is_deltasigma=True,
+            geometry=geometry,
+        )
+        assert_allclose(
+            ftDS,
+            expected["four_theta_DS"],
+            10 * reltol,
+            err_msg="4theta Shear not correct when using cluster method w/ info_mem",
+        )
+        assert_allclose(
+            cnDS,
+            expected["const_DS"],
+            10 * reltol,
+            err_msg="Constant Shear not correct when using cluster method w/ info_mem",
+        )
 
     # test basic weights functionality
     cluster.compute_galaxy_weights()
diff --git a/tests/test_galaxycluster.py b/tests/test_galaxycluster.py
index 0102cc7f0..74c124424 100644
--- a/tests/test_galaxycluster.py
+++ b/tests/test_galaxycluster.py
@@ -26,11 +26,13 @@ def test_initialization():
         "galcat": GCData(),
     }
     cl1 = clmm.GalaxyCluster(**testdict1)
+    cl2 = clmm.GalaxyCluster(**testdict1, include_quadrupole=True)
 
     assert_equal(testdict1["unique_id"], cl1.unique_id)
     assert_equal(testdict1["ra"], cl1.ra)
     assert_equal(testdict1["dec"], cl1.dec)
     assert_equal(testdict1["z"], cl1.z)
+    assert_equal(testdict1["z"], cl2.z)
     assert isinstance(cl1.galcat, GCData)
 
 
@@ -179,6 +181,13 @@ def test_integrity():  # Converge on name
     assert_equal(cl.ra, -10)
     assert_equal(cl.galcat["ra"], [-10])
 
+    # tests for phi_major
+    assert_raises(TypeError, cl.set_phi_major, phi_major=None, info_mem=None)
+    assert_raises(TypeError, cl.set_phi_major, phi_major="None", info_mem="None")
+
+    cl.set_phi_major(1.0)
+    assert_equal(cl.phi_major, 1.0)
+
 
 def test_save_load():
     """test save load"""
@@ -259,18 +268,39 @@ def test_integrity_of_lensfuncs():
     # Check metadata addition
     pzbins = np.linspace(0.0001, 5, 100)
     cluster = clmm.GalaxyCluster(unique_id="1", ra=161.3, dec=34.0, z=0.3, galcat=galcat)
+    cluster_quad = clmm.GalaxyCluster(
+        unique_id="1", ra=161.3, dec=34.0, z=0.3, galcat=galcat, include_quadrupole=True
+    )
     cluster.galcat.pzpdf_info["zbins"] = pzbins
     cluster.galcat["pzbins"] = [pzbins for i in range(len(z_src))]
     cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_src]
+    cluster_quad.galcat.pzpdf_info["zbins"] = pzbins
+    cluster_quad.galcat["pzbins"] = [pzbins for i in range(len(z_src))]
+    cluster_quad.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_src]
+
+    # check for missing quadrupole
+    assert_raises(ValueError, cluster_quad.compute_tangential_and_cross_components)
 
     for pztype in ("individual_bins", "shared_bins"):
         cluster.galcat.pzpdf_info["type"] = pztype
+        cluster_quad.galcat.pzpdf_info["type"] = pztype
 
         cluster.compute_tangential_and_cross_components(
             is_deltasigma=True, use_pdz=True, cosmo=cosmo, add=True
         )
+        cluster_quad.set_phi_major(
+            info_mem=(np.array([161.2, 161.4]), np.array([33.9, 34.1]), np.array([1.0, 1.0]))
+        )
+        cluster_quad.compute_tangential_and_cross_components(
+            is_deltasigma=True,
+            use_pdz=True,
+            cosmo=cosmo,
+            add=True,
+        )
         for comp_name in ("et", "ex"):
             assert_equal(cluster.galcat.meta[f"{comp_name}_sigmac_type"], "effective")
+        for comp_name in ("e_quad_4theta", "e_quad_const"):
+            assert_equal(cluster_quad.galcat.meta[f"{comp_name}_sigmac_type"], "effective")
 
 
 def test_integrity_of_probfuncs():
@@ -465,16 +495,35 @@ def test_plot_profiles():
             [ra_source, dec_source, shear1, shear2, z_src], names=("ra", "dec", "e1", "e2", "z")
         ),
     )
+    cluster_quad = clmm.GalaxyCluster(
+        unique_id="test",
+        ra=ra_lens,
+        dec=dec_lens,
+        z=z_lens,
+        galcat=GCData(
+            [ra_source, dec_source, shear1, shear2, z_src], names=("ra", "dec", "e1", "e2", "z")
+        ),
+        include_quadrupole=True,
+    )
     cluster.compute_tangential_and_cross_components()
+    cluster_quad.set_phi_major(
+        info_mem=(np.array([119.9, 120]), np.array([41.9, 42.1]), np.array([1.0, 1.0]))
+    )
+    cluster_quad.compute_tangential_and_cross_components()
     cluster.make_radial_profile(bin_units, bins=bins_radians, include_empty_bins=True)
+    cluster_quad.make_radial_profile(bin_units, bins=bins_radians, include_empty_bins=True)
     # missing profile name
     assert_raises(ValueError, cluster.plot_profiles, table_name="made_up_table")
+    assert_raises(ValueError, cluster_quad.plot_profiles, table_name="made_up_table")
     # missing shear component
     assert_raises(ValueError, cluster.plot_profiles, cross_component="made_up_component")
+    assert_raises(ValueError, cluster_quad.plot_profiles, quad_4theta_component="made_up_component")
     # check basic plot is working
     cluster.plot_profiles()
+    cluster_quad.plot_profiles()
     # check it passes missing a component error
     cluster.plot_profiles(cross_component_error="made_up_component")
+    cluster_quad.plot_profiles(quad_4theta_component_error="made_up_component")
 
 
 def test_coordinate_system():
diff --git a/tests/test_mockdata.py b/tests/test_mockdata.py
index 637defe75..f9e3d7599 100644
--- a/tests/test_mockdata.py
+++ b/tests/test_mockdata.py
@@ -259,13 +259,17 @@ def test_shapenoise():
 
     # Verify that the shape noise is Gaussian around 0 (for the very small shear here)
     sigma = 0.25
-    data = mock.generate_galaxy_catalog(10**12.0, 0.3, 4, cosmo, 0.8, ngals=50000, shapenoise=sigma)
+    data = mock.generate_galaxy_catalog(
+        10**12.0, 0.3, 4, cosmo, 0.8, ngals=50000, shapenoise=sigma
+    )
     # Check that there are no galaxies with |e|>1
     assert_equal(np.count_nonzero((data["e1"] > 1) | (data["e1"] < -1)), 0)
     assert_equal(np.count_nonzero((data["e2"] > 1) | (data["e2"] < -1)), 0)
     # Check that shape noise is Guassian with correct std dev
     bins = np.arange(-1, 1.1, 0.1)
-    gauss = 5000 * np.exp(-0.5 * (bins[:-1] + 0.05) ** 2 / sigma**2) / (sigma * np.sqrt(2 * np.pi))
+    gauss = (
+        5000 * np.exp(-0.5 * (bins[:-1] + 0.05) ** 2 / sigma**2) / (sigma * np.sqrt(2 * np.pi))
+    )
     assert_allclose(np.histogram(data["e1"], bins=bins)[0], gauss, atol=50, rtol=0.05)
     assert_allclose(np.histogram(data["e2"], bins=bins)[0], gauss, atol=50, rtol=0.05)
 
diff --git a/tests/test_theory.py b/tests/test_theory.py
index 99e81e057..c2f34384d 100644
--- a/tests/test_theory.py
+++ b/tests/test_theory.py
@@ -222,7 +222,7 @@ def test_compute_reduced_shear(modeling_data):
     )
 
 
-def helper_profiles(func):
+def consistency_test_density_funcs(func):
     """A helper function to repeat a set of unit tests on several functions
     that expect the same inputs.
 
@@ -275,7 +275,7 @@ def helper_profiles(func):
     )
 
 
-def test_profiles(modeling_data, profile_init):
+def test_density_profiles(modeling_data, profile_init):
     """Tests for profile functions, compute_3d_density, compute_surface_density,
     and compute_excess_surface_density"""
 
@@ -289,10 +289,10 @@ def test_profiles(modeling_data, profile_init):
         "notabackend",
         "testnotabackend",
     ]:
-        helper_profiles(theo.compute_3d_density)
-        helper_profiles(theo.compute_surface_density)
-        helper_profiles(theo.compute_mean_surface_density)
-        helper_profiles(theo.compute_excess_surface_density)
+        consistency_test_density_funcs(theo.compute_3d_density)
+        consistency_test_density_funcs(theo.compute_surface_density)
+        consistency_test_density_funcs(theo.compute_mean_surface_density)
+        consistency_test_density_funcs(theo.compute_excess_surface_density)
 
         if profile_init == "nfw":
             reltol = modeling_data["theory_reltol"]
@@ -342,97 +342,6 @@ def test_profiles(modeling_data, profile_init):
             cfg["numcosmo_profiles"]["DeltaSigma"],
             reltol,
         )
-        # Test miscentering
-        if mod.backend == "nc":
-            assert_allclose(
-                mod.eval_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                )[-40:],
-                cfg["numcosmo_profiles"]["Sigma"][-40:],
-                2.5e-2,
-            )
-            assert_allclose(
-                mod.eval_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                    mis_from_backend=True,
-                )[-40:],
-                cfg["numcosmo_profiles"]["Sigma"][-40:],
-                2.5e-2,
-            )
-            assert_allclose(
-                mod.eval_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"][-40],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                    mis_from_backend=True,
-                ),
-                cfg["numcosmo_profiles"]["Sigma"][-40],
-                2.5e-2,
-            )
-            assert_allclose(
-                mod.eval_mean_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                )[-40:],
-                (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
-                8.5e-3,
-            )
-            assert_allclose(
-                mod.eval_mean_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                    mis_from_backend=True,
-                )[-40:],
-                (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
-                8.5e-3,
-            )
-            assert_allclose(
-                mod.eval_mean_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"][-40],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                    mis_from_backend=True,
-                ),
-                (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40],
-                8.5e-3,
-            )
-            assert_allclose(
-                mod.eval_excess_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                )[-40:],
-                cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
-                2.5e-2,
-            )
-            assert_allclose(
-                mod.eval_excess_surface_density(
-                    cfg["SIGMA_PARAMS"]["r_proj"],
-                    cfg["SIGMA_PARAMS"]["z_cl"],
-                    r_mis=0.1,
-                    verbose=True,
-                    mis_from_backend=False,
-                )[-40:],
-                cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
-                2.5e-2,
-            )
-        assert_equal(theo.miscentering.integrand_surface_density_nfw(0, 0.3, 0, 0.3), 1.0 / 3.0)
-        assert_equal(
-            theo.miscentering.integrand_surface_density_hernquist(0, 0.3, 0, 0.3), 4.0 / 15.0
-        )
         if mod.backend == "ct":
             assert_raises(
                 ValueError, mod.eval_excess_surface_density, 1e-12, cfg["SIGMA_PARAMS"]["z_cl"]
@@ -469,66 +378,6 @@ def test_profiles(modeling_data, profile_init):
             reltol,
         )
 
-        # Test miscentering
-        if mod.backend == "nc":
-            assert_allclose(
-                theo.compute_surface_density(
-                    cosmo=cosmo, **cfg["SIGMA_PARAMS"], alpha_ein=alpha_ein, verbose=True, r_mis=0.1
-                )[-40:],
-                cfg["numcosmo_profiles"]["Sigma"][-40:],
-                2.5e-2,
-            )
-            assert_allclose(
-                theo.compute_surface_density(
-                    cosmo=cosmo,
-                    **cfg["SIGMA_PARAMS"],
-                    alpha_ein=alpha_ein,
-                    verbose=True,
-                    r_mis=0.1,
-                    mis_from_backend=True,
-                )[-40:],
-                cfg["numcosmo_profiles"]["Sigma"][-40:],
-                2.5e-2,
-            )
-            assert_allclose(
-                theo.compute_mean_surface_density(
-                    cosmo=cosmo, **cfg["SIGMA_PARAMS"], alpha_ein=alpha_ein, verbose=True, r_mis=0.1
-                )[-40:],
-                (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
-                8.5e-3,
-            )
-            assert_allclose(
-                theo.compute_mean_surface_density(
-                    cosmo=cosmo,
-                    **cfg["SIGMA_PARAMS"],
-                    alpha_ein=alpha_ein,
-                    verbose=True,
-                    r_mis=0.1,
-                    mis_from_backend=True,
-                )[-40:],
-                (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
-                8.5e-3,
-            )
-            assert_allclose(
-                theo.compute_excess_surface_density(
-                    cosmo=cosmo, **cfg["SIGMA_PARAMS"], alpha_ein=alpha_ein, verbose=True, r_mis=0.1
-                )[-40:],
-                cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
-                2.5e-2,
-            )
-            assert_allclose(
-                theo.compute_excess_surface_density(
-                    cosmo=cosmo,
-                    **cfg["SIGMA_PARAMS"],
-                    alpha_ein=alpha_ein,
-                    verbose=True,
-                    r_mis=0.1,
-                    mis_from_backend=True,
-                )[-40:],
-                cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
-                2.5e-2,
-            )
-
         # Test use_projected_quad
         if mod.backend == "ccl" and profile_init == "einasto":
             if hasattr(mod.hdpm, "projected_quad"):
@@ -556,45 +405,361 @@ def test_profiles(modeling_data, profile_init):
                 assert_raises(NotImplementedError, mod.set_projected_quad, True)
 
 
-def test_2halo_term(modeling_data):
+def test_miscentering(modeling_data, profile_init):
+    """Tests for profile functions, compute_3d_density, compute_surface_density,
+    and compute_excess_surface_density"""
+
+    # Validation tests
+    # NumCosmo makes different choices for constants (Msun). We make this conversion
+    # by passing the ratio of SOLAR_MASS in kg from numcosmo and CLMM
+    cfg = load_validation_config(halo_profile_model=profile_init)
+    cosmo = cfg["cosmo"]
+
+    if theo.be_nick == "nc" and modeling_data["nick"] not in [
+        "notabackend",
+        "testnotabackend",
+    ]:
+        if profile_init == "nfw":
+            reltol = modeling_data["theory_reltol"]
+        else:
+            reltol = modeling_data["theory_reltol_num"]
+
+        # Object Oriented tests
+        mod = theo.Modeling()
+        mod.set_cosmo(cosmo)
+        mod.set_halo_density_profile(halo_profile_model=cfg["SIGMA_PARAMS"]["halo_profile_model"])
+        mod.set_concentration(cfg["SIGMA_PARAMS"]["cdelta"])
+        mod.set_mass(cfg["SIGMA_PARAMS"]["mdelta"])
+
+        # Need to set the alpha value for the NC backend to the one used for the benchmarks,
+        # which is the CCL default value
+        if profile_init == "einasto":
+            alpha_ein = cfg["TEST_CASE"]["alpha_einasto"]
+            mod.set_einasto_alpha(alpha_ein)
+
+        if profile_init != "einasto":
+            alpha_ein = None
+
+        assert_equal(theo.miscentering.integrand_surface_density_nfw(0, 0.3, 0, 0.3), 1.0 / 3.0)
+        assert_equal(
+            theo.miscentering.integrand_surface_density_hernquist(0, 0.3, 0, 0.3), 4.0 / 15.0
+        )
+        if mod.backend != "nc":
+            return
+
+        # OO tests
+        assert_allclose(
+            mod.eval_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, verbose=True
+            )[-40:],
+            cfg["numcosmo_profiles"]["Sigma"][-40:],
+            2.5e-2,
+        )
+        assert_allclose(
+            mod.eval_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                r_mis=0.1,
+                verbose=True,
+                mis_from_backend=True,
+            )[-40:],
+            cfg["numcosmo_profiles"]["Sigma"][-40:],
+            2.5e-2,
+        )
+        assert_allclose(
+            mod.eval_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"][-40],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                r_mis=0.1,
+                verbose=True,
+                mis_from_backend=True,
+            ),
+            cfg["numcosmo_profiles"]["Sigma"][-40],
+            2.5e-2,
+        )
+        assert_allclose(
+            mod.eval_mean_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, verbose=True
+            )[-40:],
+            (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
+            8.5e-3,
+        )
+        assert_allclose(
+            mod.eval_mean_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                r_mis=0.1,
+                verbose=True,
+                mis_from_backend=True,
+            )[-40:],
+            (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
+            8.5e-3,
+        )
+        assert_allclose(
+            mod.eval_mean_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"][-40],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                r_mis=0.1,
+                verbose=True,
+                mis_from_backend=True,
+            ),
+            (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40],
+            8.5e-3,
+        )
+        assert_allclose(
+            mod.eval_excess_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], r_mis=0.1, verbose=True
+            )[-40:],
+            cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
+            2.5e-2,
+        )
+        assert_allclose(
+            mod.eval_excess_surface_density(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                r_mis=0.1,
+                verbose=True,
+                mis_from_backend=False,
+            )[-40:],
+            cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
+            2.5e-2,
+        )
+
+        # Function tests
+        assert_allclose(
+            theo.compute_surface_density(
+                cosmo=cosmo, **cfg["SIGMA_PARAMS"], alpha_ein=alpha_ein, verbose=True, r_mis=0.1
+            )[-40:],
+            cfg["numcosmo_profiles"]["Sigma"][-40:],
+            2.5e-2,
+        )
+        assert_allclose(
+            theo.compute_surface_density(
+                cosmo=cosmo,
+                **cfg["SIGMA_PARAMS"],
+                alpha_ein=alpha_ein,
+                verbose=True,
+                r_mis=0.1,
+                mis_from_backend=True,
+            )[-40:],
+            cfg["numcosmo_profiles"]["Sigma"][-40:],
+            2.5e-2,
+        )
+        assert_allclose(
+            theo.compute_mean_surface_density(
+                cosmo=cosmo, **cfg["SIGMA_PARAMS"], alpha_ein=alpha_ein, verbose=True, r_mis=0.1
+            )[-40:],
+            (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
+            8.5e-3,
+        )
+        assert_allclose(
+            theo.compute_mean_surface_density(
+                cosmo=cosmo,
+                **cfg["SIGMA_PARAMS"],
+                alpha_ein=alpha_ein,
+                verbose=True,
+                r_mis=0.1,
+                mis_from_backend=True,
+            )[-40:],
+            (cfg["numcosmo_profiles"]["Sigma"] + cfg["numcosmo_profiles"]["DeltaSigma"])[-40:],
+            8.5e-3,
+        )
+        assert_allclose(
+            theo.compute_excess_surface_density(
+                cosmo=cosmo, **cfg["SIGMA_PARAMS"], alpha_ein=alpha_ein, verbose=True, r_mis=0.1
+            )[-40:],
+            cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
+            2.5e-2,
+        )
+        assert_allclose(
+            theo.compute_excess_surface_density(
+                cosmo=cosmo,
+                **cfg["SIGMA_PARAMS"],
+                alpha_ein=alpha_ein,
+                verbose=True,
+                r_mis=0.1,
+                mis_from_backend=True,
+            )[-40:],
+            cfg["numcosmo_profiles"]["DeltaSigma"][-40:],
+            2.5e-2,
+        )
+
+
+def test_triaxial(modeling_data):
     cfg = load_validation_config()
     cosmo = cfg["cosmo"]
 
     # Object Oriented tests
     mod = theo.Modeling()
     mod.set_cosmo(cosmo)
+    mod.set_concentration(cfg["SIGMA_PARAMS"]["cdelta"])
+    mod.set_mass(cfg["SIGMA_PARAMS"]["mdelta"])
 
     if mod.backend not in ["ccl", "nc"]:
-        assert_raises(
-            NotImplementedError, mod.eval_surface_density_2h, 1.0, cfg["SIGMA_PARAMS"]["z_cl"]
-        )
         assert_raises(
             NotImplementedError,
-            mod.eval_excess_surface_density_2h,
+            mod.eval_excess_surface_density_triaxial,
             1.0,
             cfg["SIGMA_PARAMS"]["z_cl"],
+            0.1,
+            "mono",
+            n_grid=100,
         )
     else:
-        # Just checking that it runs and returns array of the right length
+        # Check that wrong term value is caught in OO-oriented and functional interface.
+        assert_raises(
+            ValueError,
+            theo.compute_excess_surface_density_triaxial,
+            cfg["SIGMA_PARAMS"]["r_proj"],
+            cfg["SIGMA_PARAMS"]["mdelta"],
+            cfg["SIGMA_PARAMS"]["cdelta"],
+            cfg["SIGMA_PARAMS"]["z_cl"],
+            0.1,
+            cosmo,
+            term="bleh",
+            n_grid=500,
+        )
+        assert_raises(
+            ValueError,
+            mod.eval_excess_surface_density_triaxial,
+            cfg["SIGMA_PARAMS"]["r_proj"],
+            cfg["SIGMA_PARAMS"]["z_cl"],
+            0.1,
+            "bleh",
+            n_grid=500,
+        )
+
+        # Checks that OO-oriented and functional interface give the same results
         # To be updated with proper comparison to benchmark when available
-        assert_equal(
-            len(
-                mod.eval_surface_density_2h(
-                    cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"]
-                )
+        assert_allclose(
+            theo.compute_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["mdelta"],
+                cfg["SIGMA_PARAMS"]["cdelta"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                0.1,
+                cosmo,
+                term="quad_4theta",
+                n_grid=500,
+            ),
+            mod.eval_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                0.1,
+                "quad_4theta",
+                n_grid=500,
             ),
-            len(cfg["SIGMA_PARAMS"]["r_proj"]),
+            **TOLERANCE,
         )
-        assert_equal(
-            len(
-                mod.eval_excess_surface_density_2h(
-                    cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"]
-                )
+        assert_allclose(
+            theo.compute_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["mdelta"],
+                cfg["SIGMA_PARAMS"]["cdelta"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                0.1,
+                cosmo,
+                term="quad_const",
+                n_grid=500,
             ),
-            len(cfg["SIGMA_PARAMS"]["r_proj"]),
+            mod.eval_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                0.1,
+                "quad_const",
+                n_grid=500,
+            ),
+            **TOLERANCE,
         )
+        assert_allclose(
+            theo.compute_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["mdelta"],
+                cfg["SIGMA_PARAMS"]["cdelta"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                0.1,
+                cosmo,
+                term="mono",
+                n_grid=500,
+            ),
+            mod.eval_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], 0.1, "mono", n_grid=500
+            ),
+            **TOLERANCE,
+        )
+    if mod.backend == "ccl":
+        # just to test with miscentering
+        assert_allclose(
+            theo.compute_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["mdelta"],
+                cfg["SIGMA_PARAMS"]["cdelta"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                0.1,
+                cosmo,
+                term="mono",
+                n_grid=500,
+                r_mis=1e-3,
+            ),
+            mod.eval_excess_surface_density_triaxial(
+                cfg["SIGMA_PARAMS"]["r_proj"],
+                cfg["SIGMA_PARAMS"]["z_cl"],
+                0.1,
+                "mono",
+                n_grid=500,
+                r_mis=1e-3,
+            ),
+            **TOLERANCE,
+        )
+        # just to run with verbose
+        theo.compute_excess_surface_density_triaxial(
+            cfg["SIGMA_PARAMS"]["r_proj"],
+            cfg["SIGMA_PARAMS"]["mdelta"],
+            cfg["SIGMA_PARAMS"]["cdelta"],
+            cfg["SIGMA_PARAMS"]["z_cl"],
+            0.1,
+            cosmo,
+            term="mono",
+            n_grid=20,
+            halo_profile_model="einasto",
+            verbose=True,
+        )
+        # just to run with einasto
+        theo.compute_excess_surface_density_triaxial(
+            cfg["SIGMA_PARAMS"]["r_proj"],
+            cfg["SIGMA_PARAMS"]["mdelta"],
+            cfg["SIGMA_PARAMS"]["cdelta"],
+            cfg["SIGMA_PARAMS"]["z_cl"],
+            0.1,
+            cosmo,
+            term="mono",
+            n_grid=20,
+            halo_profile_model="einasto",
+            alpha_ein=1.0,
+        )
+
+
+def test_2halo_term(modeling_data):
+    cfg = load_validation_config()
+    cosmo = cfg["cosmo"]
 
+    # Object Oriented tests
+    mod = theo.Modeling()
+    mod.set_cosmo(cosmo)
+
+    if mod.backend not in ["ccl", "nc"]:
+        assert_raises(
+            NotImplementedError, mod.eval_surface_density_2h, 1.0, cfg["SIGMA_PARAMS"]["z_cl"]
+        )
+        assert_raises(
+            NotImplementedError,
+            mod.eval_excess_surface_density_2h,
+            1.0,
+            cfg["SIGMA_PARAMS"]["z_cl"],
+        )
+    else:
         # Checks that OO-oriented and functional interface give the same results
+        # To be updated with proper comparison to benchmark when available
         assert_allclose(
             theo.compute_excess_surface_density_2h(
                 cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], cosmo
@@ -614,7 +779,7 @@ def test_2halo_term(modeling_data):
         )
 
 
-def helper_physics_functions(func, additional_kwargs={}):
+def consistency_test_shear_funcs(func, additional_kwargs={}):
     """A helper function to repeat a set of unit tests on several functions
     that expect the same inputs.
 
@@ -651,11 +816,11 @@ def helper_physics_functions(func, additional_kwargs={}):
 
 def test_shear_convergence_unittests(modeling_data, profile_init):
     """Unit and validation tests for the shear and convergence calculations"""
-    helper_physics_functions(theo.compute_tangential_shear)
-    helper_physics_functions(theo.compute_convergence)
-    helper_physics_functions(theo.compute_reduced_tangential_shear)
-    helper_physics_functions(theo.compute_magnification)
-    helper_physics_functions(theo.compute_magnification_bias, {"alpha": 1.0})
+    consistency_test_shear_funcs(theo.compute_tangential_shear)
+    consistency_test_shear_funcs(theo.compute_convergence)
+    consistency_test_shear_funcs(theo.compute_reduced_tangential_shear)
+    consistency_test_shear_funcs(theo.compute_magnification)
+    consistency_test_shear_funcs(theo.compute_magnification_bias, {"alpha": 1.0})
 
     # Validation Tests -------------------------
     # NumCosmo makes different choices for constants (Msun). We make this conversion
@@ -903,48 +1068,28 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
 
         assert_allclose(
             theo.compute_convergence(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.zeros(len(radius)),
             1.0e-10,
         )
         assert_allclose(
             theo.compute_tangential_shear(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.zeros(len(radius)),
             1.0e-10,
         )
         assert_allclose(
             theo.compute_reduced_tangential_shear(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.zeros(len(radius)),
             1.0e-10,
         )
         assert_allclose(
             theo.compute_magnification(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.ones(len(radius)),
             1.0e-10,
@@ -968,48 +1113,28 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
         z_src = [0.25, 0.1, 0.14, 0.02]
         assert_allclose(
             theo.compute_convergence(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.zeros(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
             theo.compute_tangential_shear(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.zeros(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
             theo.compute_reduced_tangential_shear(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.zeros(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
             theo.compute_magnification(
-                radius,
-                mdelta=1.0e15,
-                cdelta=4.0,
-                z_cluster=z_cluster,
-                z_src=z_src,
-                cosmo=cosmo,
+                radius, mdelta=1.0e15, cdelta=4.0, z_cluster=z_cluster, z_src=z_src, cosmo=cosmo
             ),
             np.ones(len(z_src)),
             1.0e-10,