From 09cf4d934d3e48721419374c431445c7eaaf0d85 Mon Sep 17 00:00:00 2001
From: Frank Schultz <scf175@googlemail.com>
Date: Mon, 25 Jan 2021 22:08:54 +0100
Subject: [PATCH] fix grid_lebedev() due to quadpy API changes

---
 micarray/modal/angular.py | 27 +++++++++++---
 micarray/util.py          | 76 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 99 insertions(+), 4 deletions(-)

diff --git a/micarray/modal/angular.py b/micarray/modal/angular.py
index 04afb79..7fe8922 100644
--- a/micarray/modal/angular.py
+++ b/micarray/modal/angular.py
@@ -4,7 +4,7 @@
 from .. import util
 from warnings import warn
 try:
-    import quadpy  # only for grid_lebedev()
+    import quadpy.u3._lebedev as _lebedev  # only for grid_lebedev()
 except ImportError:
     pass
 
@@ -34,6 +34,8 @@ def sht_matrix(N, azi, colat, weights=None):
     (Note: :math:`\mathbf{Y}` is interpreted as the inverse transform (or synthesis)
     matrix in examples and documentation.)
 
+    cf. eq. (1.9), (3.29) :cite:`Rafaely2019_book`
+
     Parameters
     ----------
     N : int
@@ -265,13 +267,30 @@ def available_quadrature(d):
         return matches[0]
 
     if n > 65:
+        print(n)
         raise ValueError("Maximum available Lebedev grid order is 65. "
                          "(requested: {})".format(n))
 
     # this needs https://pypi.python.org/pypi/quadpy
-    q = quadpy.sphere.Lebedev(degree=available_quadrature(2*n))
+    # currently validated with 0.16.5
+    degree = _lebedev.lebedev_003a()  # tmp call to mute pep8 warning
+    degree = available_quadrature(2 * n)
+    # the nice API handling currently returns only up to lebedev_047:
+    # q = quadpy.u3.get_good_scheme(degree=degree)
+    # thus we call the appropriate lebedev function directly in low level,
+    # not nice, but then we are save until the quadpy API will be consistent:
+    if degree == 3:
+        fcn = '_lebedev.lebedev_00' + str(degree) + 'a()'
+    elif degree < 10:
+        fcn = '_lebedev.lebedev_00'+str(degree) + '()'
+    elif degree < 100:
+        fcn = '_lebedev.lebedev_0' + str(degree) + '()'
+    else:
+        fcn = '_lebedev.lebedev_' + str(degree) + '()'
+    q = eval(fcn)
     if np.any(q.weights < 0):
         warn("Lebedev grid of order {} has negative weights.".format(n))
-    azi = q.azimuthal_polar[:, 0]
-    colat = q.azimuthal_polar[:, 1]
+    azi, colat, r = util.cart2sph(q.points[0, :],
+                                  q.points[1, :],
+                                  q.points[2, :])
     return azi, colat, 4*np.pi*q.weights
diff --git a/micarray/util.py b/micarray/util.py
index 08eb60b..2bf07b0 100644
--- a/micarray/util.py
+++ b/micarray/util.py
@@ -107,3 +107,79 @@ def db(x, power=False):
     """
     with np.errstate(divide='ignore'):
         return 10 if power else 20 * np.log10(np.abs(x))
+
+
+def sph2cart(alpha, beta, r):
+    r"""Spherical to cartesian coordinate transform.
+
+    taken from https://github.com/sfstoolbox/sfs-python commit: efdda38
+
+    .. math::
+
+        x = r \cos \alpha \sin \beta \\
+        y = r \sin \alpha \sin \beta \\
+        z = r \cos \beta
+
+    with :math:`\alpha \in [0, 2\pi), \beta \in [0, \pi], r \geq 0`
+
+    Parameters
+    ----------
+    alpha : float or array_like
+            Azimuth angle in radiants
+    beta : float or array_like
+            Colatitude angle in radiants (with 0 denoting North pole)
+    r : float or array_like
+            Radius
+
+    Returns
+    -------
+    x : float or `numpy.ndarray`
+        x-component of Cartesian coordinates
+    y : float or `numpy.ndarray`
+        y-component of Cartesian coordinates
+    z : float or `numpy.ndarray`
+        z-component of Cartesian coordinates
+
+    """
+    x = r * np.cos(alpha) * np.sin(beta)
+    y = r * np.sin(alpha) * np.sin(beta)
+    z = r * np.cos(beta)
+    return x, y, z
+
+
+def cart2sph(x, y, z):
+    r"""Cartesian to spherical coordinate transform.
+
+    taken from https://github.com/sfstoolbox/sfs-python commit: efdda38
+
+    .. math::
+
+        \alpha = \arctan \left( \frac{y}{x} \right) \\
+        \beta = \arccos \left( \frac{z}{r} \right) \\
+        r = \sqrt{x^2 + y^2 + z^2}
+
+    with :math:`\alpha \in [-pi, pi], \beta \in [0, \pi], r \geq 0`
+
+    Parameters
+    ----------
+    x : float or array_like
+        x-component of Cartesian coordinates
+    y : float or array_like
+        y-component of Cartesian coordinates
+    z : float or array_like
+        z-component of Cartesian coordinates
+
+    Returns
+    -------
+    alpha : float or `numpy.ndarray`
+            Azimuth angle in radiants
+    beta : float or `numpy.ndarray`
+            Colatitude angle in radiants (with 0 denoting North pole)
+    r : float or `numpy.ndarray`
+            Radius
+
+    """
+    r = np.sqrt(x**2 + y**2 + z**2)
+    alpha = np.arctan2(y, x)
+    beta = np.arccos(z / r)
+    return alpha, beta, r