From 3dd3d64c44d94821e6f3457a23c61e994d01dccb Mon Sep 17 00:00:00 2001
From: Radhakrishnan Srinivasan
 <sriniva@x1403c0s6b0n1h0.chn.perlmutter.nersc.gov>
Date: Wed, 17 Jul 2024 05:08:38 -0700
Subject: [PATCH 1/8] Minor comments and cosmetic changres

---
 clmm/theory/func_layer.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/clmm/theory/func_layer.py b/clmm/theory/func_layer.py
index d69e629a0..699d462e0 100644
--- a/clmm/theory/func_layer.py
+++ b/clmm/theory/func_layer.py
@@ -1308,7 +1308,7 @@ def compute_delta_sigma_4theta_triaxiality(
     verbose=False,
     validate_input=True,
 ):
-    r"""Compute the "4-theta"-quadrupole value
+    r"""Compute the "4-theta" component of quadrupole shear as given in Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366)
 
     Parameters
     ----------
@@ -1351,7 +1351,7 @@ def compute_delta_sigma_4theta_triaxiality(
     Returns
     -------
     ds4theta: array
-        Delta sigma 4-theta value at a position for the elliptical halo specified
+        Delta sigma 4-theta component value at a position for the elliptical halo specified
     """
     ### DEFINING INTEGRALS:
     r_arr = np.linspace(0.01, 3*np.max(r_source), sample_N)
@@ -1390,7 +1390,7 @@ def compute_delta_sigma_const_triaxiality(
     verbose=False,
     validate_input=True,
 ):
-    r"""Compute the excess surface density lensing profile for "const"-quadrupole value
+    r"""Compute the "const" component of quadrupole shear as given in Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366)
 
     Parameters
     ----------
@@ -1474,7 +1474,8 @@ def compute_delta_sigma_excess_triaxiality(
     verbose=False,
     validate_input=True,
 ):
-    r"""Compute the excess surface density lensing profile for the monopole component along with second order expansion term (e**2)
+    r"""Compute the excess surface density lensing profile for the monopole component along 
+    with second order expansion term (e**2) as given in Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366)
 
     Parameters
     ----------

From 2f3997927dbbda0b4f694b16c6016bc51f30a978 Mon Sep 17 00:00:00 2001
From: Radhakrishnan Srinivasan
 <sriniva@x1403c0s6b0n1h0.chn.perlmutter.nersc.gov>
Date: Wed, 17 Jul 2024 05:53:46 -0700
Subject: [PATCH 2/8] Data operation-functions for triaxiality added to dataops

---
 clmm/__init__.py         |   2 +-
 clmm/dataops/__init__.py | 187 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 188 insertions(+), 1 deletion(-)

diff --git a/clmm/__init__.py b/clmm/__init__.py
index ee4a5e9c3..7d98303e4 100644
--- a/clmm/__init__.py
+++ b/clmm/__init__.py
@@ -2,7 +2,7 @@
 from .gcdata import GCData
 from .galaxycluster import GalaxyCluster
 from .clusterensemble import ClusterEnsemble
-from .dataops import compute_tangential_and_cross_components, make_radial_profile
+from .dataops import compute_tangential_and_cross_components, make_radial_profile, make_binned_estimators_triaxiality
 from .utils import compute_radial_averages, make_bins, convert_units
 from .theory import (
     compute_reduced_shear_from_convergence,
diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index a93cba312..8a0e16f69 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -2,6 +2,7 @@
 import warnings
 import numpy as np
 import scipy
+from scipy.stats import binned_statistic
 from astropy.coordinates import SkyCoord
 from astropy import units as u
 from ..gcdata import GCData
@@ -708,3 +709,189 @@ def make_stacked_radial_profile(angsep, weights, components):
         np.average(component, axis=0, weights=weights) for component in components
     ]
     return staked_angsep, stacked_components
+
+def measure_Delta_Sigma_const_triaxiality(w, gamma1, Sigma_crit) :
+    """Measure Delta sigma const quadrupole value for individual galaxies.
+
+    Parameters
+    ----------
+    w: float or array
+        Weight computed for each galaxy by the descrition as given in Shin et al. (https://doi.org/10.1093/mnras/stx3366)
+    gamma1: float or array
+        Shear component gamma 1
+    Sigma_crit: float
+        Critical surface mass density in M_{sun}/Mpc^{2}
+
+    Returns
+    -------
+    dsconst_data: The measured "4theta" quadrupole component from the input shear measurements for each galaxy
+    """
+    dsconst_data = w * Sigma_crit * gamma1 / w
+    return dsconst_data
+
+def measure_Delta_Sigma_4theta_triaxiality(w1, w2, gamma1, gamma2, theta, Sigma_crit) :
+    """Measure Delta sigma 4theta quadrupole value for individual galaxies.
+
+    Parameters
+    ----------
+    w1: float or array
+        Weight computed for each galaxy by the descrition as given in Shin et al. (https://doi.org/10.1093/mnras/stx3366)
+    w2: float or array
+        Weight computed for each galaxy by the descrition as given in Shin et al. (https://doi.org/10.1093/mnras/stx3366)
+    gamma1: float or array
+        Shear component gamma 1
+    gamma2: float or array
+        Shear component gamma 2
+    theta: float or array
+        Position Angle of galaxy measured counter-clockwise from the gamma1 direction
+    Sigma_crit: float
+        Critical surface mass density in M_{sun}/Mpc^{2}
+
+    Returns
+    -------
+    ds4theta_data: The measured "const" quadrupole component from the input shear measurements for each galaxy
+    """
+    ds4theta_data = Sigma_crit * (w1*gamma1/np.cos(4*theta) + w2*gamma2/np.sin(4*theta)) / (w1 + w2)
+    return ds4theta_data
+
+def measure_weights_triaxiality(Sigma_crit, theta, Sigma_shape=0.0001, Sigma_meas=0) :
+    """Measure weight values for quadrupole measurement of individual galaxies. 
+    Using Equations 35, 32, 33 from Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366)
+
+    Parameters
+    ----------
+    Sigma_crit: float
+        Critical surface mass density in M_{sun}/Mpc^{2}
+    theta: float or array
+        Position Angle of galaxy measured counter-clockwise from the gamma1 direction
+    Sigma_shape: Float, optional
+        Scatter in the shape measurement from shape noise
+    Sigma_meas: Float, optional
+        Measurement error for shapes
+
+    Returns
+    -------
+    w, w1, w2: float or array
+        Weights for each galaxy
+    """
+	
+    w  = 1 / (Sigma_crit**2 * (Sigma_shape**2 + Sigma_meas**2))
+    w1 = np.cos(4*theta)**2 * w
+    w2 = np.sin(4*theta)**2 * w
+    return w, w1, w2
+
+def make_binned_estimators_triaxiality(gamma1, gamma2, x_arcsec, y_arcsec, z_cluster, z_source, cosmo, monopole_bins=15,
+                                      quadrupole_bins=15, r_min=0.3, r_max=2.5):
+    """
+    Makes radially binned profiles for Monopole, Quadrupole (4theta and const) as described in Shin et al. 2018
+    Parameters
+    ----------
+    gamma1: float or array
+        Shear component gamma 1
+    gamma2: float or array
+        Shear component gamma 2
+    x_arcsec: float or array
+        X position in arcsec when cluster center is at 0.0 and x axis is assumed to be along the major axis of the cluster
+    y_arcsec: float or array
+        Y position in arcsec when cluster center is at 0.0
+    z_cluster: float
+        Redshift of lens cluster
+    cosmo: clmm.cosmology.Cosmology object
+        CLMM Cosmology object 
+    monopole_bins: float, optional
+        Number of bins for monopole lensing profile
+    quadrupole_bins: float, optional
+        Number of bins for quadrupole lensing profiles (4theta, const)
+    r_min: Float, optional
+        Minimum radial distance to measure monopole and quadrupole shear profiles in Mpc (Default is 0.3 Mpc)
+    r_max: Float, optional
+        Maximum radial distance to measure monopole and quadrupole shear profiles in Mpc (Default is 2.5 Mpc)
+
+    Returns
+    -------
+    ds_mono: float or array
+        Binned monopole lensing profile
+    ds_mono_err: float or array
+        Uncertainty for the binned monopole lensing shear profile
+    r_mono: float or array
+        Bin centers for the monopole lensing shear profile
+    ds_4theta: float or array
+        Binned binned quadrupole 4-theta lensing shear profile
+    ds_4theta_err: float or array
+        Uncertainty for the binned quadrupole 4-theta lensing shear profile
+    ds_const: float or array
+        Binned binned quadrupole const lensing shear profile
+    ds_const_err: float or array
+        Uncertainty for the binned quadrupole const lensing shear profile
+    r_quad: float or array
+        Bin centers for the quadrupole lensing shear profiles (4theta and const)
+    """
+    
+    sigma_crit = cosmo.eval_sigma_crit(z_cluster,z_source)
+    r = np.sqrt((x_arcsec**2 + y_arcsec**2)) #In arcsecs
+    theta = np.arctan2(y_arcsec, x_arcsec)
+    r_mpc = r*cosmo.eval_da(z_cluster) * np.pi/180.0 * 1/3600 #In Mpc
+    
+    w, w1, w2 = measure_weights_triaxiality(sigma_crit, theta)
+    DS4theta = measure_Delta_Sigma_4theta_triaxiality(w1, w2, gamma1, gamma2, theta, sigma_crit)
+    DSconst = measure_Delta_Sigma_const_triaxiality(w, gamma1, sigma_crit)
+    
+    #Binning for Quadrupole measurements DS4theta and DSconst
+    bins=quadrupole_bins
+    r_min = r_min
+    r_max = r_max
+    
+    bin_edges = np.logspace(np.log10(r_min), np.log10(r_max), bins)
+    N_i = []
+    for i in np.arange(bins-1):
+        N_i.append(len(r_mpc[(r_mpc > bin_edges[i]) & (r_mpc < bin_edges[i+1])]))
+    N_i=np.array(N_i)
+
+
+    result = binned_statistic(r_mpc, gamma1, statistic='mean', bins=bin_edges)
+    gamma1_i = result.statistic
+    res = binned_statistic(r_mpc, gamma2, statistic='mean', bins=bin_edges)
+    gamma2_i = res.statistic
+    res = binned_statistic(r_mpc, DS4theta, statistic='mean', bins=bin_edges)
+    DS4theta_i_err = binned_statistic(r_mpc, DS4theta, statistic='std', bins=bin_edges).statistic/np.sqrt(N_i)
+    DS4theta_i = res.statistic
+    res = binned_statistic(r_mpc, DSconst, statistic='mean', bins=bin_edges)
+    DSconst_i = res.statistic
+    DSconst_i_err = binned_statistic(r_mpc, DSconst, statistic='std', bins=bin_edges).statistic/np.sqrt(N_i)
+    r_i = bin_edges
+    
+    #Binning for Monopole Measurements:
+    bins_mono=monopole_bins
+    bin_edges = np.logspace(np.log10(r_min), np.log10(r_max), bins_mono)
+    N_i = []
+    for i in np.arange(bins_mono-1):
+        N_i.append(len(r_mpc[(r_mpc > bin_edges[i]) & (r_mpc < bin_edges[i+1])]))
+    N_i=np.array(N_i)
+
+    r_mono = bin_edges
+    res = binned_statistic(r_mpc, -gamma1*np.cos(2*theta)-gamma2*np.sin(2*theta), statistic='mean', bins=bin_edges)
+    gammat_mono = res.statistic
+    ds_mono_err = binned_statistic(r_mpc, -gamma1*np.cos(2*theta)-gamma2*np.sin(2*theta), statistic='std',
+                                   bins=bin_edges).statistic/np.sqrt(N_i)*sigma_crit
+    ds_mono = gammat_mono*sigma_crit
+        
+    # SAFEGUARD AGAINST BINS WITH NANs and 0.0s
+    
+    ind = np.invert(np.isnan(ds_mono) | np.isnan(ds_mono_err))
+    ds_mono = ds_mono[ind]
+    ds_mono_err = ds_mono_err[ind]
+    r_mono = np.sqrt(r_mono[:-1]*r_mono[1:])[ind]
+    ind = (ds_mono!= 0.0) & (ds_mono_err!= 0.0)
+    ds_mono = ds_mono[ind]
+    ds_mono_err = np.abs(ds_mono_err[ind]) 
+    r_mono = r_mono[ind] 
+    
+    ind = np.invert(np.isnan(DS4theta_i) | np.isnan(DS4theta_i_err) | np.isnan(DSconst_i) | np.isnan(DSconst_i_err))
+    ds_4theta = DS4theta_i[ind]
+    ds_4theta_err = np.abs(DS4theta_i_err[ind])
+    ds_const = DSconst_i[ind] 
+    ds_const_err = np.abs(DSconst_i_err[ind])
+    r_quad = np.sqrt(r_i[:-1]*r_i[1:])[ind]
+    
+    return ds_mono,ds_mono_err,r_mono,ds_4theta,ds_4theta_err,ds_const,ds_const_err,r_quad
+

From 46067b883d95329946bdbc676920a3e22caea392 Mon Sep 17 00:00:00 2001
From: tae-h-shin <babohahaman@gmail.com>
Date: Wed, 17 Jul 2024 09:13:24 -0400
Subject: [PATCH 3/8] edited galaxycluster.py to incorporate quadrupole lensing

---
 clmm/galaxycluster.py | 340 ++++++++++++++++++++++++++++++++----------
 1 file changed, 263 insertions(+), 77 deletions(-)

diff --git a/clmm/galaxycluster.py b/clmm/galaxycluster.py
index b4db5446c..7564a82fc 100644
--- a/clmm/galaxycluster.py
+++ b/clmm/galaxycluster.py
@@ -33,31 +33,57 @@ class GalaxyCluster:
         Declination of galaxy cluster center (in degrees)
     z : float
         Redshift of galaxy cluster center
+    phi_major : float
+        Direction of the cluster major axis (in radian).
+        Users may provide `ra_mem`, `dec_mem` and `weight_mem` instead of this parameter. 
+        Only needed when `is_quadrupole==True`.
+    ra_mem : array_like
+        RAs of member galaxies (in degrees).
+        Only needed when `is_quadrupole==True`.
+    dec_mem : array_like
+        DECs of member galaxies (in degrees).
+        Only needed when `is_quadrupole==True`.
+    weight_mem : array_like
+        weights of member galaxies (to calculate major axis).
+        Only needed when `is_quadrupole==True`.
     galcat : GCData
         Table of background galaxy data containing at least galaxy_id, ra, dec, e1, e2, z
     validate_input: bool
         Validade each input argument
+    is_quadrupole: bool
+        If quadrupole WL be calculated.
     """
 
-    def __init__(self, *args, validate_input=True, **kwargs):
+    def __init__(self, *args, validate_input=True, is_quadrupole=False, **kwargs):
         self.unique_id = None
         self.ra = None
         self.dec = None
         self.z = None
         self.galcat = None
+        self.phi_major = None
+        self.ra_mem = None
+        self.dec_mem = None
+        self.weight_mem = None
         self.validate_input = validate_input
+        self.is_quadrupole = is_quadrupole
         if len(args) > 0 or len(kwargs) > 0:
             self._add_values(*args, **kwargs)
             self._check_types()
             self.set_ra_lower(ra_low=0)
 
-    def _add_values(self, unique_id: str, ra: float, dec: float, z: float, galcat: GCData):
+    def _add_values(self, unique_id: str, ra: float, dec: float, z: float, galcat: GCData, 
+                    phi_major=None, ra_mem=None, dec_mem=None, weight_mem=None):
         """Add values for all attributes"""
         self.unique_id = unique_id
         self.ra = ra
         self.dec = dec
         self.z = z
         self.galcat = galcat
+        if self.is_quadrupole:
+            self.phi_major = phi_major
+            self.ra_mem = ra_mem
+            self.dec_mem = dec_mem
+            self.weight_mem = weight_mem
 
     def _check_types(self):
         """Check types of all attributes"""
@@ -70,6 +96,17 @@ def _check_types(self):
         self.ra = float(self.ra)
         self.dec = float(self.dec)
         self.z = float(self.z)
+        if self.is_quadrupole:
+            if self.phi_major is not None:
+                validate_argument(vars(self), "phi_major", float)
+                self.phi_major = float(self.phi_major)
+            else:
+                validate_argument(vars(self), "ra_mem", "float_array")
+                validate_argument(vars(self), "dec_mem", "float_array")
+                validate_argument(vars(self), "weight_mem", "float_array")
+                self.ra_mem =  np.array(self.ra_mem, dtype=float)
+                self.dec_mem = np.array(self.dec_mem, dtype=float)
+                self.weight_mem = np.array(self.weight_mem, dtype=float)
 
     def save(self, filename, **kwargs):
         """Saves GalaxyCluster object to filename using Pickle"""
@@ -208,6 +245,8 @@ def compute_tangential_and_cross_components(
         shape_component2="e2",
         tan_component="et",
         cross_component="ex",
+        quad_4theta_component="e4theta",
+        quad_const_component="econst",
         geometry="curve",
         is_deltasigma=False,
         use_pdz=False,
@@ -225,6 +264,13 @@ def compute_tangential_and_cross_components(
         shear2: `galcat` shape_component2
         geometry: `input` geometry
         is_deltasigma: `input` is_deltasigma
+        is_quadrupole: `input` is_quadrupole
+        [if is_quadrupole:]
+            phi_major: `cluster` major axis direction (in radian with respect to +x)
+            (OR)
+            ra_mem: `cluster` RAs of member galaxies for calculating major axis of a given cluster
+            dec_mem: `cluster` DECs of member galaxies for calculating major axis of a given cluster
+            weight_mem: `cluster` weights to be applied to member galaxies for calculating major axis
 
         Parameters
         ----------
@@ -242,6 +288,14 @@ 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: `e4theta`
+        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: `econst`
         geometry: str, optional
             Sky geometry to compute angular separation.
             Options are curve (uses astropy) or flat.
@@ -257,10 +311,16 @@ def compute_tangential_and_cross_components(
         -------
         angsep: array_like
             Angular separation between lens and each source galaxy in radians
-        tangential_component: array_like
-            Tangential shear (or assimilated quantity) for each source galaxy
-        cross_component: array_like
-            Cross shear (or assimilated quantity) for each source galaxy
+        [if not self.is_quadrupole:]
+            tangential_component: array_like
+                Tangential shear (or assimilated quantity) for each source galaxy
+            cross_component: array_like
+                Cross shear (or assimilated quantity) for each source galaxy
+        [if self.is_quadrupole:]
+            quad_4theta_component: array_like
+                4theta quadrupole shear (or assimilated quantity) for each source galaxy
+            quad_const_component: array_like
+                constatnt quadrupole shear (or assimilated quantity) for each source galaxy
         """
         # Check is all the required data is available
         col_dict = {
@@ -275,23 +335,47 @@ 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(
-            is_deltasigma=is_deltasigma,
-            ra_lens=self.ra,
-            dec_lens=self.dec,
-            geometry=geometry,
-            validate_input=self.validate_input,
-            **cols,
-        )
-        if add:
-            self.galcat["theta"] = angsep
-            self.galcat[tan_component] = tangential_comp
-            self.galcat[cross_component] = cross_comp
-            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
+        if self.is_quadrupole:
+            angsep, four_theta_comp, const_comp = compute_tangential_and_cross_components(
+                is_deltasigma=is_deltasigma,
+                ra_lens=self.ra,
+                dec_lens=self.dec,
+                geometry=geometry,
+                validate_input=self.validate_input,
+                is_quadupole=self.is_quarupole,
+                phi_major=self.phi_major,
+                ra_mem=self.ra_mem,
+                dec_mem=self.dec_mem,
+                weight_mem=self.weight_mem,
+                **cols,
+            )
+            if add:
+                self.galcat["theta"] = angsep
+                self.galcat[quad_4theta_component] = four_theta_comp
+                self.galcat[quad_const_component] = const_comp
+                if is_deltasigma:
+                    sigmac_type = "effective" if use_pdz else "standard"
+                    self.galcat.meta[f"{quad_4theta_component}_sigmac_type"] = sigmac_type
+                    self.galcat.meta[f"{quad_const_component}_sigmac_type"] = sigmac_type
+            return angsep, four_theta_comp, const_comp
+        else:
+            angsep, tangential_comp, cross_comp = compute_tangential_and_cross_components(
+                is_deltasigma=is_deltasigma,
+                ra_lens=self.ra,
+                dec_lens=self.dec,
+                geometry=geometry,
+                validate_input=self.validate_input,
+                **cols,
+            )
+            if add:
+                self.galcat["theta"] = angsep
+                self.galcat[tan_component] = tangential_comp
+                self.galcat[cross_component] = cross_comp
+                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
 
     def compute_background_probability(
         self, use_pdz=False, add=True, p_background_name="p_background"
@@ -468,10 +552,16 @@ def make_radial_profile(
         cosmo=None,
         tan_component_in="et",
         cross_component_in="ex",
+        quad_4theta_component_in="e4theta",
+        quad_const_component_in="econst",
         tan_component_out="gt",
         cross_component_out="gx",
+        quad_4theta_component_out="g4theta",
+        quad_const_component_out="gconst",
         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,
@@ -488,6 +578,8 @@ def make_radial_profile(
 
         Calls `clmm.dataops.make_radial_profile` with the following arguments:
         components: `galcat` components (tan_component_in, cross_component_in, z)
+                    OR
+                    (quad_4theta_component_in, quad_const_component_in, z) IF is_quadrupole
         angsep: `galcat` theta
         angsep_units: 'radians'
         bin_units: `input` bin_units
@@ -521,18 +613,36 @@ 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: 'e4theta'
+        quad_const_component_in: string, optional
+            Name of the constant quadrupole component column in `galcat` to be binned.
+            Default: 'econst'
         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: 'g4theta'
+        quad_const_component_out: string, optional
+            Name of the constant quadrupole component binned profile column to be added in profile table.
+            Default: 'gconst'
         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
@@ -562,39 +672,73 @@ 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")
-        ):
-            raise TypeError(
-                "Shear or ellipticity information is missing. Galaxy catalog must have tangential"
-                "and cross shears (gt, gx) or ellipticities (et, ex). "
-                "Run compute_tangential_and_cross_components first."
-            )
+        if self.is_quadrupole:
+            if not all(
+                t_ in self.galcat.columns for t_ in (quad_4theta_component_in, quad_const_component_in, "theta")
+            ):
+                raise TypeError(
+                    "Shear or ellipticity information is missing. Galaxy catalog must have 4theta"
+                    "and constant quadrupole shears (g4theta, gconst) or ellipticities (e4theta, econst). "
+                    "Run compute_tangential_and_cross_components first."
+                )
+        else:
+            if not all(
+                t_ in self.galcat.columns for t_ in (tan_component_in, cross_component_in, "theta")
+            ):
+                raise TypeError(
+                    "Shear or ellipticity information is missing. Galaxy catalog must have tangential"
+                    "and cross shears (gt, gx) or ellipticities (et, ex). "
+                    "Run compute_tangential_and_cross_components first."
+                )
         if "z" not in self.galcat.columns:
             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")],
-            angsep=self.galcat["theta"],
-            angsep_units="radians",
-            bin_units=bin_units,
-            bins=bins,
-            error_model=error_model,
-            include_empty_bins=include_empty_bins,
-            return_binnumber=True,
-            cosmo=cosmo,
-            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)
-            ],
-            weights=self.galcat[weights_in].data if use_weights else None,
-        )
-        # Reaname table columns
-        for i, name in enumerate([tan_component_out, cross_component_out, "z"]):
-            profile_table.rename_column(f"p_{i}", name)
-            profile_table.rename_column(f"p_{i}_err", f"{name}_err")
+        if self.is_quadrupole:
+            profile_table, binnumber = make_radial_profile(
+                [self.galcat[n].data for n in (quad_4theta_component_in, quad_const_component_in, "z")],
+                angsep=self.galcat["theta"],
+                angsep_units="radians",
+                bin_units=bin_units,
+                bins=bins,
+                error_model=error_model,
+                include_empty_bins=include_empty_bins,
+                return_binnumber=True,
+                cosmo=cosmo,
+                z_lens=self.z,
+                validate_input=self.validate_input,
+                components_error=[
+                    None if n is None else self.galcat[n].data
+                    for n in (quad_4theta_component_in_err, quad_const_component_in_err, None)
+                ],
+                weights=self.galcat[weights_in].data if use_weights else None,
+            )
+            # Reaname table columns
+            for i, name in enumerate([quad_4theta_component_out, quad_const_component_out, "z"]):
+                profile_table.rename_column(f"p_{i}", name)
+                profile_table.rename_column(f"p_{i}_err", f"{name}_err")
+        else:
+            profile_table, binnumber = make_radial_profile(
+                [self.galcat[n].data for n in (tan_component_in, cross_component_in, "z")],
+                angsep=self.galcat["theta"],
+                angsep_units="radians",
+                bin_units=bin_units,
+                bins=bins,
+                error_model=error_model,
+                include_empty_bins=include_empty_bins,
+                return_binnumber=True,
+                cosmo=cosmo,
+                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)
+                ],
+                weights=self.galcat[weights_in].data if use_weights else None,
+            )
+            # Reaname table columns
+            for i, name in enumerate([tan_component_out, cross_component_out, "z"]):
+                profile_table.rename_column(f"p_{i}", name)
+                profile_table.rename_column(f"p_{i}_err", f"{name}_err")
         # Reaname weights columns
         profile_table.rename_column("weights_sum", weights_out)
         # add galaxy IDs
@@ -626,6 +770,10 @@ def plot_profiles(
         tangential_component_error="gt_err",
         cross_component="gx",
         cross_component_error="gx_err",
+        quad_4theta_component="g4theta",
+        quad_4theta_component_error="g4theta_err",
+        quad_const_component="gconst",
+        quad_const_component_error="gconst_err",
         table_name="profile",
         xscale="linear",
         yscale="linear",
@@ -646,6 +794,18 @@ 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: 'g4theta'
+        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: 'g4theta_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:
@@ -663,30 +823,56 @@ 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,
-        )
+        if self.is_quadrupole:
+            for col in (quad_4theta_component, quad_const_component):
+                if col not in profile.columns:
+                    raise ValueError(f"Column for plotting '{col}' does not exist.")
+            for col in (quad_4theta_component_error, quad_const_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[quad_4theta_component],
+                tangential_component_error=(
+                    profile[quad_4theta_component_error]
+                    if quad_4theta_component_error in profile.columns
+                    else None
+                ),
+                cross_component=profile[quad_const_component],
+                cross_component_error=(
+                    profile[quad_const_component_error] if quad_const_component_error in profile.columns else None
+                ),
+                xscale=xscale,
+                yscale=yscale,
+                tangential_component_label=quad_4theta_component,
+                cross_component_label=quad_const_component,
+            )
+        else:
+            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,
+            )
 
     def set_ra_lower(self, ra_low=0):
         """

From 8cef454d0eb17e4aa3a7cf2f9905a63b5faf385c Mon Sep 17 00:00:00 2001
From: tae-h-shin <babohahaman@gmail.com>
Date: Wed, 17 Jul 2024 09:18:07 -0400
Subject: [PATCH 4/8] edited galaxycluster.py to incorporate quadrupole lensing

---
 examples/demo_mock_ensemble.ipynb | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/examples/demo_mock_ensemble.ipynb b/examples/demo_mock_ensemble.ipynb
index 84354b111..a14eb088f 100644
--- a/examples/demo_mock_ensemble.ipynb
+++ b/examples/demo_mock_ensemble.ipynb
@@ -618,9 +618,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "wrk",
+   "display_name": "Python [conda env:tf-gpu]",
    "language": "python",
-   "name": "wrk"
+   "name": "conda-env-tf-gpu-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -632,7 +632,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.9.18"
   }
  },
  "nbformat": 4,

From 8525ba7cd715359c2069d6da51ffd43ed067bb4b Mon Sep 17 00:00:00 2001
From: tae-h-shin <babohahaman@gmail.com>
Date: Wed, 17 Jul 2024 10:08:12 -0400
Subject: [PATCH 5/8] changed dataops/__init__.py for quadrupole

---
 clmm/dataops/__init__.py | 53 +++++++++++++++++++++-------------------
 1 file changed, 28 insertions(+), 25 deletions(-)

diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index 8a0e16f69..e42482cea 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -33,7 +33,7 @@ def compute_tangential_and_cross_components(
     geometry="curve",
     is_deltasigma=False,
     sigma_c=None,
-    is_quadrupole=False,
+    include_quadrupole=False,
     phi_major=None,
     ra_mem=None,
     dec_mem=None,
@@ -118,21 +118,21 @@ 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.
-    is_quadrupole: bool
+    include_quadrupole: bool
         If `True`, the quadrupole shear components (g_4theta, g_const; Shin+2018) are calculated
         instead of g_t and g_x
     phi_major : float, optional
         the direction of the major axis of the input cluster in the unit of radian. 
-        only needed when `is_quadupole` is `True`.
+        only needed when `include_quadrupole` is `True`.
         Users could choose to provide ra_mem, dec_mem and weight_mem instead of this quantity.
     ra_mem : array, optional
         right ascentions of the member galaxies of the input cluster,
         to calculate the direction of the major axis of the cluster.
-        Only needed when `is_quadrupole` is `True` and `phi_major` is not provided.
+        Only needed when `include_quadrupole` is `True` and `phi_major` is not provided.
     dec_mem : array, optional
         declinations of the member galaxies of the input cluster,
         to calculate the direction of the major axis of the cluster.
-        Only needed when `is_quadrupole` is `True` and `phi_major` is not provided.
+        Only needed when `include_quadrupole` is `True` and `phi_major` is not provided.
     weight_mem : array, optional
         weights for the member galaxies to be applied when calcualting 
         the major axis of the input cluster out of ra_mem and dec_mem. 
@@ -144,16 +144,15 @@ def compute_tangential_and_cross_components(
     -------
     angsep: array_like
         Angular separation between lens and each source galaxy in radians
-    IF is_quadrupole:
+    tangential_component: array_like
+        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
-    ELSE:
-        tangential_component: array_like
-            Tangential shear (or assimilated quantity) for each source galaxy
-        cross_component: array_like
-            Cross shear (or assimilated quantity) for each source galaxy
     """
     # pylint: disable-msg=too-many-locals
     # Note: we make these quantities to be np.array so that a name is not passed from astropy
@@ -167,7 +166,7 @@ def compute_tangential_and_cross_components(
         validate_argument(locals(), "shear2", "float_array")
         validate_argument(locals(), "geometry", str)
         validate_argument(locals(), "is_deltasigma", bool)
-        validate_argument(locals(), "is_quadrupole", bool)
+        validate_argument(locals(), "include_quadrupole", bool)
         validate_argument(locals(), "sigma_c", "float_array", none_ok=True)
         ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency(
             [ra_source, dec_source, shear1, shear2],
@@ -175,6 +174,12 @@ def compute_tangential_and_cross_components(
             prefix="Tangential- and Cross- shape components sources",
         )
         _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c)
+        if phi_major is not None:
+            validate_argument(locals(), "phi_major", float)
+        else:
+            validate_argument(locals(), "ra_mem", "float_array")
+            validate_argument(locals(), "dec_mem", "float_array")
+            validate_argument(locals(), "weight_mem", "float_array")
     elif np.iterable(ra_source):
         ra_source_, dec_source_, shear1_, shear2_ = (
             np.array(col) for col in [ra_source, dec_source, shear1, shear2]
@@ -193,7 +198,15 @@ def compute_tangential_and_cross_components(
         angsep, phi = _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_, dec_source_)
     else:
         raise NotImplementedError(f"Sky geometry {geometry} is not currently supported")
-    if is_quadupole:
+    # Compute the tangential and cross shears
+    tangential_comp = _compute_tangential_shear(shear1_, shear2_, phi)
+    cross_comp = _compute_cross_shear(shear1_, shear2_, phi)
+    # 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
+    if include_quadrupole:
         if phi_major is not None:
             phi_major_ = phi_major
         else:
@@ -207,20 +220,10 @@ def compute_tangential_and_cross_components(
         const_comp = rotated_shear1
         # 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)
             four_theta_comp *= _sigma_c_arr
             const_comp *= _sigma_c_arr
-    else:
-        # Compute the tangential and cross shears
-        tangential_comp = _compute_tangential_shear(shear1_, shear2_, phi)
-        cross_comp = _compute_cross_shear(shear1_, shear2_, phi)
-        # 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
-    if is_quadrupole:
-        return angsep, four_theta_comp, const_comp
+    if include_quadrupole:
+        return angsep, tangential_comp, cross_comp, four_theta_comp, const_comp
     else:
         return angsep, tangential_comp, cross_comp
 

From 40358a7434d068f5b60a2b7a8fbd8b8c233a7216 Mon Sep 17 00:00:00 2001
From: tae-h-shin <babohahaman@gmail.com>
Date: Wed, 17 Jul 2024 10:30:33 -0400
Subject: [PATCH 6/8] changed galaxycluser.py for quadrupole

---
 clmm/dataops/__init__.py |   4 +-
 clmm/galaxycluster.py    | 109 +++++++++++++++++++++++++--------------
 2 files changed, 72 insertions(+), 41 deletions(-)

diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index e42482cea..602dcde9a 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -206,6 +206,7 @@ def compute_tangential_and_cross_components(
         _sigma_c_arr = np.array(sigma_c)
         tangential_comp *= _sigma_c_arr
         cross_comp *= _sigma_c_arr
+    
     if include_quadrupole:
         if phi_major is not None:
             phi_major_ = phi_major
@@ -221,12 +222,13 @@ def compute_tangential_and_cross_components(
         # If the is_deltasigma flag is True, multiply the results by Sigma_crit.
         if sigma_c is not None:
             four_theta_comp *= _sigma_c_arr
-            const_comp *= _sigma_c_arr
+            const_comp *= _sigma_c_arr    
     if include_quadrupole:
         return angsep, tangential_comp, cross_comp, four_theta_comp, const_comp
     else:
         return angsep, tangential_comp, cross_comp
 
+
 def compute_background_probability(
     z_lens, z_src=None, use_pdz=False, pzpdf=None, pzbins=None, validate_input=True
 ):
diff --git a/clmm/galaxycluster.py b/clmm/galaxycluster.py
index 7564a82fc..35a109c65 100644
--- a/clmm/galaxycluster.py
+++ b/clmm/galaxycluster.py
@@ -36,25 +36,25 @@ class GalaxyCluster:
     phi_major : float
         Direction of the cluster major axis (in radian).
         Users may provide `ra_mem`, `dec_mem` and `weight_mem` instead of this parameter. 
-        Only needed when `is_quadrupole==True`.
+        Only needed when `include_quadrupole==True`.
     ra_mem : array_like
         RAs of member galaxies (in degrees).
-        Only needed when `is_quadrupole==True`.
+        Only needed when `include_quadrupole==True`.
     dec_mem : array_like
         DECs of member galaxies (in degrees).
-        Only needed when `is_quadrupole==True`.
+        Only needed when `include_quadrupole==True`.
     weight_mem : array_like
         weights of member galaxies (to calculate major axis).
-        Only needed when `is_quadrupole==True`.
+        Only needed when `include_quadrupole==True`.
     galcat : GCData
         Table of background galaxy data containing at least galaxy_id, ra, dec, e1, e2, z
     validate_input: bool
         Validade each input argument
-    is_quadrupole: bool
+    include_quadrupole: bool
         If quadrupole WL be calculated.
     """
 
-    def __init__(self, *args, validate_input=True, is_quadrupole=False, **kwargs):
+    def __init__(self, *args, validate_input=True, include_quadrupole=False, **kwargs):
         self.unique_id = None
         self.ra = None
         self.dec = None
@@ -65,7 +65,7 @@ def __init__(self, *args, validate_input=True, is_quadrupole=False, **kwargs):
         self.dec_mem = None
         self.weight_mem = None
         self.validate_input = validate_input
-        self.is_quadrupole = is_quadrupole
+        self.include_quadrupole = include_quadrupole
         if len(args) > 0 or len(kwargs) > 0:
             self._add_values(*args, **kwargs)
             self._check_types()
@@ -79,7 +79,7 @@ def _add_values(self, unique_id: str, ra: float, dec: float, z: float, galcat: G
         self.dec = dec
         self.z = z
         self.galcat = galcat
-        if self.is_quadrupole:
+        if self.include_quadrupole:
             self.phi_major = phi_major
             self.ra_mem = ra_mem
             self.dec_mem = dec_mem
@@ -96,7 +96,7 @@ def _check_types(self):
         self.ra = float(self.ra)
         self.dec = float(self.dec)
         self.z = float(self.z)
-        if self.is_quadrupole:
+        if self.include_quadrupole:
             if self.phi_major is not None:
                 validate_argument(vars(self), "phi_major", float)
                 self.phi_major = float(self.phi_major)
@@ -264,8 +264,8 @@ def compute_tangential_and_cross_components(
         shear2: `galcat` shape_component2
         geometry: `input` geometry
         is_deltasigma: `input` is_deltasigma
-        is_quadrupole: `input` is_quadrupole
-        [if is_quadrupole:]
+        include_quadrupole: `input` include_quadrupole
+        [if include_quadrupole:]
             phi_major: `cluster` major axis direction (in radian with respect to +x)
             (OR)
             ra_mem: `cluster` RAs of member galaxies for calculating major axis of a given cluster
@@ -311,12 +311,11 @@ def compute_tangential_and_cross_components(
         -------
         angsep: array_like
             Angular separation between lens and each source galaxy in radians
-        [if not self.is_quadrupole:]
-            tangential_component: array_like
-                Tangential shear (or assimilated quantity) for each source galaxy
-            cross_component: array_like
-                Cross shear (or assimilated quantity) for each source galaxy
-        [if self.is_quadrupole:]
+        tangential_component: array_like
+            Tangential shear (or assimilated quantity) for each source galaxy
+        cross_component: array_like
+            Cross shear (or assimilated quantity) for each source galaxy
+        [if self.include_quadrupole:]
             quad_4theta_component: array_like
                 4theta quadrupole shear (or assimilated quantity) for each source galaxy
             quad_const_component: array_like
@@ -335,14 +334,14 @@ def compute_tangential_and_cross_components(
         cols = self._get_input_galdata(col_dict)
 
         # compute shears
-        if self.is_quadrupole:
-            angsep, four_theta_comp, const_comp = compute_tangential_and_cross_components(
+        if self.include_quadrupole:
+            angsep, tangential_comp, cross_comp, four_theta_comp, const_comp = compute_tangential_and_cross_components(
                 is_deltasigma=is_deltasigma,
                 ra_lens=self.ra,
                 dec_lens=self.dec,
                 geometry=geometry,
                 validate_input=self.validate_input,
-                is_quadupole=self.is_quarupole,
+                include_quadrupole=self.include_quadrupole,
                 phi_major=self.phi_major,
                 ra_mem=self.ra_mem,
                 dec_mem=self.dec_mem,
@@ -351,13 +350,17 @@ def compute_tangential_and_cross_components(
             )
             if add:
                 self.galcat["theta"] = angsep
+                self.galcat[tan_component] = tangential_comp
+                self.galcat[cross_component] = cross_comp
                 self.galcat[quad_4theta_component] = four_theta_comp
                 self.galcat[quad_const_component] = const_comp
                 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
                     self.galcat.meta[f"{quad_4theta_component}_sigmac_type"] = sigmac_type
                     self.galcat.meta[f"{quad_const_component}_sigmac_type"] = sigmac_type
-            return angsep, four_theta_comp, const_comp
+            return angsep, tangential_comp, cross_comp, four_theta_comp, const_comp
         else:
             angsep, tangential_comp, cross_comp = compute_tangential_and_cross_components(
                 is_deltasigma=is_deltasigma,
@@ -579,7 +582,7 @@ def make_radial_profile(
         Calls `clmm.dataops.make_radial_profile` with the following arguments:
         components: `galcat` components (tan_component_in, cross_component_in, z)
                     OR
-                    (quad_4theta_component_in, quad_const_component_in, z) IF is_quadrupole
+                    (quad_4theta_component_in, quad_const_component_in, z) IF include_quadrupole
         angsep: `galcat` theta
         angsep_units: 'radians'
         bin_units: `input` bin_units
@@ -671,8 +674,15 @@ def make_radial_profile(
         """
         # Too many local variables (19/15)
         # pylint: disable=R0914
-
-        if self.is_quadrupole:
+        if not all(
+            t_ in self.galcat.columns for t_ in (tan_component_in, cross_component_in, "theta")
+        ):
+            raise TypeError(
+                "Shear or ellipticity information is missing. Galaxy catalog must have tangential"
+                "and cross shears (gt, gx) or ellipticities (et, ex). "
+                "Run compute_tangential_and_cross_components first."
+            )
+        if self.include_quadrupole:
             if not all(
                 t_ in self.galcat.columns for t_ in (quad_4theta_component_in, quad_const_component_in, "theta")
             ):
@@ -681,21 +691,16 @@ def make_radial_profile(
                     "and constant quadrupole shears (g4theta, gconst) or ellipticities (e4theta, econst). "
                     "Run compute_tangential_and_cross_components first."
                 )
-        else:
-            if not all(
-                t_ in self.galcat.columns for t_ in (tan_component_in, cross_component_in, "theta")
-            ):
-                raise TypeError(
-                    "Shear or ellipticity information is missing. Galaxy catalog must have tangential"
-                    "and cross shears (gt, gx) or ellipticities (et, ex). "
-                    "Run compute_tangential_and_cross_components first."
-                )
         if "z" not in self.galcat.columns:
             raise TypeError("Missing galaxy redshifts!")
         # Compute the binned averages and associated errors
-        if self.is_quadrupole:
+        if self.include_quadrupole:
             profile_table, binnumber = make_radial_profile(
-                [self.galcat[n].data for n in (quad_4theta_component_in, quad_const_component_in, "z")],
+                [self.galcat[n].data for n in 
+                    (tan_component_in, cross_component_in, 
+                        quad_4theta_component_in, quad_const_component_in, "z"
+                    )
+                ],
                 angsep=self.galcat["theta"],
                 angsep_units="radians",
                 bin_units=bin_units,
@@ -708,12 +713,16 @@ def make_radial_profile(
                 validate_input=self.validate_input,
                 components_error=[
                     None if n is None else self.galcat[n].data
-                    for n in (quad_4theta_component_in_err, quad_const_component_in_err, None)
+                    for n in (tan_component_in_err, cross_component_in_err, 
+                        quad_4theta_component_in_err, quad_const_component_in_err, None
+                    )
                 ],
                 weights=self.galcat[weights_in].data if use_weights else None,
             )
             # Reaname table columns
-            for i, name in enumerate([quad_4theta_component_out, quad_const_component_out, "z"]):
+            for i, name in enumerate([tan_component_out, cross_component_out,
+                quad_4theta_component_out, quad_const_component_out, "z"]
+            ):
                 profile_table.rename_column(f"p_{i}", name)
                 profile_table.rename_column(f"p_{i}_err", f"{name}_err")
         else:
@@ -823,14 +832,33 @@ 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)
-        if self.is_quadrupole:
-            for col in (quad_4theta_component, quad_const_component):
+        if self.include_quadrupole:
+            for col in (tangential_component, cross_component, quad_4theta_component, quad_const_component):
                 if col not in profile.columns:
                     raise ValueError(f"Column for plotting '{col}' does not exist.")
-            for col in (quad_4theta_component_error, quad_const_component_error):
+            for col in (tangential_component_error, cross_component_error, 
+                    quad_4theta_component_error, quad_const_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,
+            ), plot_profiles(
                 rbins=profile["radius"],
                 r_units=profile.meta["bin_units"],
                 tangential_component=profile[quad_4theta_component],
@@ -848,6 +876,7 @@ def plot_profiles(
                 tangential_component_label=quad_4theta_component,
                 cross_component_label=quad_const_component,
             )
+
         else:
             for col in (tangential_component, cross_component):
                 if col not in profile.columns:

From ad70d4ed3f64d058150e6ccb59cc1bf1770888ee Mon Sep 17 00:00:00 2001
From: Radhakrishnan Srinivasan
 <sriniva@x1403c0s3b0n0h0.chn.perlmutter.nersc.gov>
Date: Thu, 18 Jul 2024 01:28:52 -0700
Subject: [PATCH 7/8] Edited a minor typo in code in dataops/__init__.py

---
 clmm/dataops/__init__.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index 8a0e16f69..fdfef169c 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -123,7 +123,7 @@ def compute_tangential_and_cross_components(
         instead of g_t and g_x
     phi_major : float, optional
         the direction of the major axis of the input cluster in the unit of radian. 
-        only needed when `is_quadupole` is `True`.
+        only needed when `is_quadrupole` is `True`.
         Users could choose to provide ra_mem, dec_mem and weight_mem instead of this quantity.
     ra_mem : array, optional
         right ascentions of the member galaxies of the input cluster,
@@ -193,7 +193,7 @@ def compute_tangential_and_cross_components(
         angsep, phi = _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_, dec_source_)
     else:
         raise NotImplementedError(f"Sky geometry {geometry} is not currently supported")
-    if is_quadupole:
+    if is_quadrupole:
         if phi_major is not None:
             phi_major_ = phi_major
         else:
@@ -690,7 +690,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
@@ -726,7 +726,7 @@ def measure_Delta_Sigma_const_triaxiality(w, gamma1, Sigma_crit) :
     -------
     dsconst_data: The measured "4theta" quadrupole component from the input shear measurements for each galaxy
     """
-    dsconst_data = w * Sigma_crit * gamma1 / w
+    dsconst_data = Sigma_crit * gamma1 
     return dsconst_data
 
 def measure_Delta_Sigma_4theta_triaxiality(w1, w2, gamma1, gamma2, theta, Sigma_crit) :

From 2a196f7586641e5219fc3f8f6c43fb4856969ef9 Mon Sep 17 00:00:00 2001
From: Radhakrishnan Srinivasan
 <sriniva@x1403c0s3b0n0h0.chn.perlmutter.nersc.gov>
Date: Thu, 18 Jul 2024 01:38:11 -0700
Subject: [PATCH 8/8] Edited a small typo in dataops/__init__.py

---
 clmm/dataops/__init__.py | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index 875f2e7f6..cd04cd60b 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -123,11 +123,7 @@ def compute_tangential_and_cross_components(
         instead of g_t and g_x
     phi_major : float, optional
         the direction of the major axis of the input cluster in the unit of radian. 
-<<<<<<< HEAD
-        only needed when `is_quadrupole` is `True`.
-=======
         only needed when `include_quadrupole` is `True`.
->>>>>>> 40358a7434d068f5b60a2b7a8fbd8b8c233a7216
         Users could choose to provide ra_mem, dec_mem and weight_mem instead of this quantity.
     ra_mem : array, optional
         right ascentions of the member galaxies of the input cluster,
@@ -202,9 +198,6 @@ def compute_tangential_and_cross_components(
         angsep, phi = _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_, dec_source_)
     else:
         raise NotImplementedError(f"Sky geometry {geometry} is not currently supported")
-<<<<<<< HEAD
-    if is_quadrupole:
-=======
     # Compute the tangential and cross shears
     tangential_comp = _compute_tangential_shear(shear1_, shear2_, phi)
     cross_comp = _compute_cross_shear(shear1_, shear2_, phi)
@@ -215,7 +208,6 @@ def compute_tangential_and_cross_components(
         cross_comp *= _sigma_c_arr
     
     if include_quadrupole:
->>>>>>> 40358a7434d068f5b60a2b7a8fbd8b8c233a7216
         if phi_major is not None:
             phi_major_ = phi_major
         else: