Repository: skypyproject/skypy
Pull request #407 (2021-01-22 17:30:07)

random projection of triaxial ellipsoid
## Description

Factors out the random projection of the ellipsoid from `skypy.galaxy.ellipticity.ryden04` into a random sampler `skypy.utils.random.triaxial_axis_ratio`.

This might be reused for other ellipticity sampling for galaxies (@jucordero) or halos and clusters (@sutieng).

## Checklist
- [x] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/master/CONTRIBUTING.md)
- [x] Write unit tests
- [x] Write documentation strings
- [x] Assign someone from your working team to review this pull request
- [x] Assign someone from the infrastructure team to review this pull request



Patch:
diff --git a/skypy/galaxy/ellipticity.py b/skypy/galaxy/ellipticity.py
index 1cf6997..af2fe8f 100644
--- a/skypy/galaxy/ellipticity.py
+++ b/skypy/galaxy/ellipticity.py
@@ -5,6 +5,7 @@ This module provides facilities to sample the ellipticities of galaxies.
 
 import numpy as np
 from scipy import stats
+from ..utils import random
 
 
 __all__ = [
@@ -120,19 +121,8 @@ def ryden04(mu_gamma, sigma_gamma, mu, sigma, size=None):
     if size == () and not np.isscalar(gam):  # pragma: no cover
         gam, eps = gam.item(), eps.item()
 
-    # draw random viewing angle (theta, phi)
-    cos2_theta = np.random.uniform(low=-1., high=1., size=size)**2
-    cos2_phi = np.cos(np.random.uniform(low=0., high=2*np.pi, size=size))**2
-    sin2_theta = 1 - cos2_theta
-    sin2_phi = 1 - cos2_phi
-
-    # compute Binney's ABC -- eq.s (13)-(15)
-    A = (1 - eps*(2-eps)*sin2_phi)*cos2_theta + gam**2*sin2_theta
-    B = (2*eps*(2-eps))**2*cos2_theta*sin2_phi*cos2_phi
-    C = 1 - eps*(2-eps)*cos2_phi
-
-    # compute axis ratio q -- eq. (12)
-    q = np.sqrt((A+C-np.sqrt((A-C)**2+B))/(A+C+np.sqrt((A-C)**2+B)))
+    # random projection of random triaxial ellipsoid
+    q = random.triaxial_axis_ratio(1-eps, gam)
 
     # return the ellipticity
     return (1-q)/(1+q)
diff --git a/skypy/utils/random.py b/skypy/utils/random.py
index 0686a31..30500a0 100644
--- a/skypy/utils/random.py
+++ b/skypy/utils/random.py
@@ -11,6 +11,7 @@ Utility Functions
    :toctree: ../api/
 
    schechter
+   triaxial_axis_ratio
 
 """
 
@@ -73,3 +74,61 @@ def schechter(alpha, x_min, x_max, resolution=100, size=None, scale=1.):
     lnx_sample = np.interp(u, cdf, lnx)
 
     return np.exp(lnx_sample) * scale
+
+
+def triaxial_axis_ratio(zeta, xi, size=None):
+    r'''axis ratio of a randomly projected triaxial ellipsoid
+
+    Given the two axis ratios `1 >= zeta >= xi` of a randomly oriented triaxial
+    ellipsoid, computes the axis ratio `q` of the projection.
+
+    Parameters
+    ----------
+    zeta : array_like
+        Axis ratio of intermediate and major axis.
+    xi : array_like
+        Axis ratio of minor and major axis.
+    size : tuple of int or None
+        Size of the random draw. If `None` is given, size is inferred from
+        other inputs.
+
+    Returns
+    -------
+    q : array_like
+        Axis ratio of the randomly projected ellipsoid.
+
+    Notes
+    -----
+    See equations (11) and (12) in [1]_ for details.
+
+    References
+    ----------
+    .. [1] Binney J., 1985, MNRAS, 212, 767. doi:10.1093/mnras/212.4.767
+    '''
+
+    # get size from inputs if not explicitly provided
+    if size is None:
+        size = np.broadcast(zeta, xi).shape
+
+    # draw random viewing angle (theta, phi)
+    cos2_theta = np.random.uniform(low=-1., high=1., size=size)
+    cos2_theta *= cos2_theta
+    sin2_theta = 1 - cos2_theta
+    cos2_phi = np.cos(np.random.uniform(low=0., high=2*np.pi, size=size))
+    cos2_phi *= cos2_phi
+    sin2_phi = 1 - cos2_phi
+
+    # transform arrays to quantities that are used in eq. (11)
+    z2m1 = np.square(zeta)
+    z2m1 -= 1
+    x2 = np.square(xi)
+
+    # eq. (11) multiplied by xi^2 zeta^2
+    A = (1 + z2m1*sin2_phi)*cos2_theta + x2*sin2_theta
+    B2 = 4*z2m1**2*cos2_theta*sin2_phi*cos2_phi
+    C = 1 + z2m1*cos2_phi
+
+    # eq. (12)
+    q = np.sqrt((A+C-np.sqrt((A-C)**2+B2))/(A+C+np.sqrt((A-C)**2+B2)))
+
+    return q


Test patch:
diff --git a/skypy/utils/tests/test_random.py b/skypy/utils/tests/test_random.py
index 4f8a706..534200c 100644
--- a/skypy/utils/tests/test_random.py
+++ b/skypy/utils/tests/test_random.py
@@ -56,3 +56,37 @@ def test_schechter_gamma():
 
     _, p = kstest(x, 'gamma', args=(alpha+1, 0, scale))
     assert p > 0.01
+
+
+def test_triaxial_axis_ratio():
+
+    from skypy.utils.random import triaxial_axis_ratio
+
+    # sample a single axis ratio
+    q = triaxial_axis_ratio(0.8, 0.4)
+    assert np.isscalar(q)
+
+    # sample many axis ratios
+    q = triaxial_axis_ratio(0.8, 0.4, size=1000)
+    assert np.shape(q) == (1000,)
+
+    # sample with explicit shape
+    q = triaxial_axis_ratio(0.8, 0.4, size=(10, 10))
+    assert np.shape(q) == (10, 10)
+
+    # sample with implicit size
+    q1 = triaxial_axis_ratio([0.8, 0.9], 0.4)
+    q2 = triaxial_axis_ratio(0.8, [0.4, 0.5])
+    assert np.shape(q1) == np.shape(q2) == (2,)
+
+    # sample with broadcasting rule
+    q = triaxial_axis_ratio([[0.6, 0.7], [0.8, 0.9]], [0.4, 0.5])
+    assert np.shape(q) == (2, 2)
+
+    # sample with random parameters and check that projection is
+    # between largest and smallest possible value
+    zeta, xi = np.sort(np.random.uniform(0, 1, size=(2, 1000)), axis=0)
+    qmin = np.min([zeta, xi, xi/zeta], axis=0)
+    qmax = np.max([zeta, xi, xi/zeta], axis=0)
+    q = triaxial_axis_ratio(zeta, xi)
+    assert np.all((qmax >= q) & (q >= qmin))