feat: Add coordinate_system parameter to compute_tangential_and_cross_components and _compute_lensing_angles functions

clmm/__init__.py

@@ -26,4 +26,4 @@
 )
 from . import support
 
-__version__ = "1.12.3"
+__version__ = "1.12.4"

clmm/dataops/__init__.py

@@ -1,4 +1,5 @@
 """Functions to compute polar/azimuthal averages in radial bins"""
+
 import warnings
 import numpy as np
 import scipy
@@ -14,6 +15,7 @@
     _validate_ra,
     _validate_dec,
     _validate_is_deltasigma_sigma_c,
+    _validate_coordinate_system,
 )
 from ..redshift import (
     _integ_pzfuncs,
@@ -29,6 +31,7 @@ def compute_tangential_and_cross_components(
     dec_source,
     shear1,
     shear2,
+    coordinate_system="euclidean",
     geometry="curve",
     is_deltasigma=False,
     sigma_c=None,
@@ -95,6 +98,10 @@ def compute_tangential_and_cross_components(
         The measured shear (or reduced shear or ellipticity) of the source galaxies
     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.
+        Default is 'euclidean'.
     geometry: str, optional
         Sky geometry to compute angular separation.
         Options are curve (uses astropy) or flat.
@@ -128,6 +135,7 @@ def compute_tangential_and_cross_components(
         validate_argument(locals(), "shear2", "float_array")
         validate_argument(locals(), "geometry", str)
         validate_argument(locals(), "sigma_c", "float_array", none_ok=True)
+        _validate_coordinate_system(locals(), "coordinate_system", str)
         ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency(
             [ra_source, dec_source, shear1, shear2],
             names=("Ra", "Dec", "Shear1", "Shear2"),
@@ -147,9 +155,13 @@ def compute_tangential_and_cross_components(
         )
     # Compute the lensing angles
     if geometry == "flat":
-        angsep, phi = _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_, dec_source_)
+        angsep, phi = _compute_lensing_angles_flatsky(
+            ra_lens, dec_lens, ra_source_, dec_source_, coordinate_system=coordinate_system
+        )
     elif geometry == "curve":
-        angsep, phi = _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_, dec_source_)
+        angsep, phi = _compute_lensing_angles_astropy(
+            ra_lens, dec_lens, ra_source_, dec_source_, coordinate_system=coordinate_system
+        )
     else:
         raise NotImplementedError(f"Sky geometry {geometry} is not currently supported")
     # Compute the tangential and cross shears
@@ -329,7 +341,9 @@ def compute_galaxy_weights(
     return w_ls
 
 
-def _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_list, dec_source_list):
+def _compute_lensing_angles_flatsky(
+    ra_lens, dec_lens, ra_source_list, dec_source_list, coordinate_system="euclidean"
+):
     r"""Compute the angular separation between the lens and the source and the azimuthal
     angle from the lens to the source in radians.
 
@@ -353,6 +367,10 @@ def _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_list, dec_sourc
         Right ascensions of each source galaxy in degrees
     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'.
 
     Returns
     -------
@@ -376,12 +394,16 @@ def _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_list, dec_sourc
         phi[angsep == 0.0] = 0.0
     else:
         phi = 0.0 if angsep == 0.0 else phi
+    if coordinate_system == "celestial":
+        phi = np.pi - phi
     if np.any(angsep > np.pi / 180.0):
         warnings.warn("Using the flat-sky approximation with separations >1 deg may be inaccurate")
     return angsep, phi
 
 
-def _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_list, dec_source_list):
+def _compute_lensing_angles_astropy(
+    ra_lens, dec_lens, ra_source_list, dec_source_list, coordinate_system="euclidean"
+):
     r"""Compute the angular separation between the lens and the source and the azimuthal
     angle from the lens to the source in radians.
 
@@ -395,6 +417,10 @@ def _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_list, dec_sourc
         Right ascensions of each source galaxy in degrees
     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'.
 
     Returns
     -------
@@ -414,6 +440,8 @@ def _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_list, dec_sourc
     else:
         phi -= 2 * np.pi if phi > np.pi else 0
         phi = 0 if angsep == 0 else phi
+    if coordinate_system == "celestial":
+        phi = np.pi - phi
     return angsep, phi
 
 

clmm/galaxycluster.py

@@ -1,6 +1,7 @@
 """@file galaxycluster.py
 The GalaxyCluster class
 """
+
 import pickle
 import warnings
 from .gcdata import GCData
@@ -17,6 +18,7 @@
     _validate_ra,
     _validate_dec,
     _draw_random_points_from_tab_distribution,
+    _validate_coordinate_system,
 )
 
 
@@ -35,6 +37,10 @@ class GalaxyCluster:
         Redshift of galaxy cluster center
     galcat : GCData
         Table of background galaxy data containing at least galaxy_id, ra, dec, e1, e2, z
+    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'.
     validate_input: bool
         Validade each input argument
     """
@@ -51,13 +57,22 @@ def __init__(self, *args, validate_input=True, **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,
+        coordinate_system: str = "euclidean",
+    ):
         """Add values for all attributes"""
         self.unique_id = unique_id
         self.ra = ra
         self.dec = dec
         self.z = z
         self.galcat = galcat
+        self.coordinate_system = coordinate_system
 
     def _check_types(self):
         """Check types of all attributes"""
@@ -66,6 +81,7 @@ def _check_types(self):
         _validate_dec(vars(self), "dec", False)
         validate_argument(vars(self), "z", (float, str), argmin=0, eqmin=True)
         validate_argument(vars(self), "galcat", GCData)
+        _validate_coordinate_system(vars(self), "coordinate_system", str)
         self.unique_id = str(self.unique_id)
         self.ra = float(self.ra)
         self.dec = float(self.dec)
@@ -281,6 +297,7 @@ def compute_tangential_and_cross_components(
             dec_lens=self.dec,
             geometry=geometry,
             validate_input=self.validate_input,
+            coordinate_system=self.coordinate_system,
             **cols,
         )
         if add:

clmm/support/mock_data.py

@@ -1,4 +1,5 @@
 """Functions to generate mock source galaxy distributions to demo lensing code"""
+
 import warnings
 import numpy as np
 from astropy import units as u
@@ -12,6 +13,7 @@
     validate_argument,
     _draw_random_points_from_distribution,
     gaussian,
+    _validate_coordinate_system,
 )
 from ..redshift import distributions as zdist
 
@@ -38,6 +40,7 @@ def generate_galaxy_catalog(
     ngal_density=None,
     pz_bins=101,
     pzpdf_type="shared_bins",
+    coordinate_system="euclidean",
     validate_input=True,
 ):
     r"""Generates a mock dataset of sheared background galaxies.
@@ -156,6 +159,11 @@ def generate_galaxy_catalog(
         The number density of galaxies (in galaxies per square arcminute, from z=0 to z=infty).
         The number of galaxies to be drawn will then depend on the redshift distribution and
         user-defined redshift range.  If specified, the ngals argument will be ignored.
+    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'.
+
     validate_input: bool
         Validade each input argument
 
@@ -208,6 +216,7 @@ def generate_galaxy_catalog(
         validate_argument(locals(), "ngals", float, none_ok=True)
         validate_argument(locals(), "ngal_density", float, none_ok=True)
         validate_argument(locals(), "pz_bins", (int, "array"))
+        _validate_coordinate_system(locals(), "coordinate_system", str)
 
     if zsrc_min is None:
         zsrc_min = cluster_z + 0.1
@@ -231,6 +240,7 @@ def generate_galaxy_catalog(
         "pz_bins": pz_bins,
         "field_size": field_size,
         "pzpdf_type": pzpdf_type,
+        "coordinate_system": coordinate_system,
     }
 
     if ngals is None and ngal_density is None:
@@ -361,6 +371,7 @@ def _generate_galaxy_catalog(
     photoz_sigma_unscaled=None,
     pz_bins=101,
     pzpdf_type="shared_bins",
+    coordinate_system="euclidean",
     field_size=None,
 ):
     """A private function that skips the sanity checks on derived properties. This
@@ -453,6 +464,9 @@ def _generate_galaxy_catalog(
             cols += ["pzbins"]
         cols += ["pzpdf"]
 
+    if coordinate_system == "celestial":
+        galaxy_catalog["e2"] *= -1  # flip e2 to match the celestial coordinate system
+
     return galaxy_catalog[cols]
 
 

clmm/utils/__init__.py

@@ -40,6 +40,7 @@
     _validate_ra,
     _validate_dec,
     _validate_is_deltasigma_sigma_c,
+    _validate_coordinate_system,
 )
 
 from .units import (

clmm/utils/validation.py

@@ -1,4 +1,5 @@
 """General utility functions that are used in multiple modules"""
+
 import numpy as np
 from ..constants import Constants as const
 
@@ -224,3 +225,23 @@ def _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c):
         raise TypeError("sigma_c (=None) must be provided when is_deltasigma=True")
     if not is_deltasigma and sigma_c is not None:
         raise TypeError(f"sigma_c (={sigma_c}) must be None when is_deltasigma=False")
+
+
+def _validate_coordinate_system(loc, coordinate_system, valid_type):
+    r"""Validate the coordinate system.
+
+    Parameters
+    ----------
+    loc: dict
+        Dictionary with all input arguments. Should be locals().
+    coordinate_system: str
+        Coordinate system of the ellipticity components. Must be either 'celestial' or 'euclidean'.
+    valid_type: str, type
+        Valid types for argument, options are object types, list/tuple of types, or:
+
+            * 'int_array' - interger, interger array
+            * 'float_array' - float, float array
+    """
+    validate_argument(loc, coordinate_system, valid_type)
+    if loc[coordinate_system] not in ["celestial", "euclidean"]:
+        raise ValueError(f"{coordinate_system} must be 'celestial' or 'euclidean'.")

examples/demo_dataops_functionality.ipynb

@@ -106,7 +106,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Define toy cluster parameters for mock data generation"
+    "Define toy cluster parameters for mock data generation. Notice that here we set the coordinate system for the ellipticity components as 'euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details."
    ]
   },
   {
@@ -139,6 +139,7 @@
     "    ngals=ngals,\n",
     "    cluster_ra=cluster_ra,\n",
     "    cluster_dec=cluster_dec,\n",
+    "    coordinate_system=\"euclidean\",\n",
     ")"
    ]
   },
@@ -155,7 +156,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cl = GalaxyCluster(cluster_id, cluster_ra, cluster_dec, cluster_z, noisy_data_z)"
+    "cl = GalaxyCluster(\n",
+    "    cluster_id, cluster_ra, cluster_dec, cluster_z, noisy_data_z, coordinate_system=\"euclidean\"\n",
+    ")"
    ]
   },
   {
@@ -252,6 +255,13 @@
     "### 3.1 Compute angular separation, cross and tangential shear for each source galaxy"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To keep things consistent, you also need to set the coordinate system to the same option."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -265,6 +275,7 @@
     "    dec_source=cl.galcat[\"dec\"],\n",
     "    shear1=cl.galcat[\"e1\"],\n",
     "    shear2=cl.galcat[\"e2\"],\n",
+    "    coordinate_system=\"euclidean\",\n",
     ")"
    ]
   },
@@ -550,7 +561,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cl.make_radial_profile(\"Mpc\", cosmo=cosmo, gal_ids_in_bins=True);"
+    "cl.make_radial_profile(\"Mpc\", cosmo=cosmo, gal_ids_in_bins=True)"
    ]
   },
   {
@@ -704,6 +715,7 @@
     "                       tan_component_out='DeltaSigma_tan', cross_component_out='DeltaSigma_cross',\n",
     "                       table_name='DeltaSigma_profile').pprint(max_width=-1)\n",
     "\"\"\"\n",
+    "\n",
     "cl.make_radial_profile(\n",
     "    \"Mpc\",\n",
     "    cosmo=cosmo,\n",
@@ -713,7 +725,7 @@
     "    cross_component_out=\"DeltaSigma_cross\",\n",
     "    table_name=\"DeltaSigma_profile\",\n",
     "    use_weights=True,\n",
-    ");"
+    ")"
    ]
   },
   {
@@ -796,7 +808,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.10"
+   "version": "3.12.4"
   }
  },
  "nbformat": 4,