Python scipy.special.sph_harm() Examples

def sph_harm(m, n, az, el):
    """Compute spherical harmonics

    Parameters
    ----------
    m : (int)
        Order of the spherical harmonic. abs(m) <= n

    n : (int)
        Degree of the harmonic, sometimes called l. n >= 0

    az: (float)
        Azimuthal (longitudinal) coordinate [0, 2pi], also called Theta.

    el : (float)
        Elevation (colatitudinal) coordinate [0, pi], also called Phi.

    Returns
    -------
    y_mn : (complex float)
        Complex spherical harmonic of order m and degree n,
        sampled at theta = az, phi = el
    """

    return scy.sph_harm(m, n, az, el) 

def process_coords_for_computations(self, coords_for_computations, s, t):
        """
        """
        if self._teffext:
            return coords_for_computations

        x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1))
        theta = np.arccos(z/r)
        phi = np.arctan2(y, x)

        xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
        xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
        xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)

        new_coords = np.zeros(coords_for_computations.shape)
        new_coords[:,0] = coords_for_computations[:,0] + xi_r * np.sin(theta) * np.cos(phi)
        new_coords[:,1] = coords_for_computations[:,1] + xi_r * np.sin(theta) * np.sin(phi)
        new_coords[:,2] = coords_for_computations[:,2] + xi_r * np.cos(theta)

        return new_coords 

def sph_harm_all(nMax, az, el):
    """Compute all spherical harmonic coefficients up to degree nMax.

    Parameters
    ----------
    nMax : (int)
        Maximum degree of coefficients to be returned. n >= 0

    az: (float), array_like
        Azimuthal (longitudinal) coordinate [0, 2pi], also called Theta.

    el : (float), array_like
        Elevation (colatitudinal) coordinate [0, pi], also called Phi.

    Returns
    -------
    y_mn : (complex float), array_like
        Complex spherical harmonics of degrees n [0 ... nMax] and all corresponding
        orders m [-n ... n], sampled at [az, el]. dim1 corresponds to az/el pairs,
        dim2 to oder/degree (m, n) pairs like 0/0, -1/1, 0/1, 1/1, -2/2, -1/2 ...
    """
    m, n = mnArrays(nMax)
    mA, azA = _np.meshgrid(m, az)
    nA, elA = _np.meshgrid(n, el)
    return sph_harm(mA, nA, azA, elA) 

def test_sph_harm():
    # Tests derived from tables in
    # http://en.wikipedia.org/wiki/Table_of_spherical_harmonics
    sh = special.sph_harm
    pi = np.pi
    exp = np.exp
    sqrt = np.sqrt
    sin = np.sin
    cos = np.cos
    yield (assert_array_almost_equal, sh(0,0,0,0),
           0.5/sqrt(pi))
    yield (assert_array_almost_equal, sh(-2,2,0.,pi/4),
           0.25*sqrt(15./(2.*pi)) *
           (sin(pi/4))**2.)
    yield (assert_array_almost_equal, sh(-2,2,0.,pi/2),
           0.25*sqrt(15./(2.*pi)))
    yield (assert_array_almost_equal, sh(2,2,pi,pi/2),
           0.25*sqrt(15/(2.*pi)) *
           exp(0+2.*pi*1j)*sin(pi/2.)**2.)
    yield (assert_array_almost_equal, sh(2,4,pi/4.,pi/3.),
           (3./8.)*sqrt(5./(2.*pi)) *
           exp(0+2.*pi/4.*1j) *
           sin(pi/3.)**2. *
           (7.*cos(pi/3.)**2.-1))
    yield (assert_array_almost_equal, sh(4,4,pi/8.,pi/6.),
           (3./16.)*sqrt(35./(2.*pi)) *
           exp(0+4.*pi/8.*1j)*sin(pi/6.)**4.) 

def Ylm(l, m, theta, phi):
    return sph_harm(m, l, phi, theta) 

def process_coords_for_observations(self, coords_for_computations, coords_for_observations, s, t):
        """
        Displacement equations:

          xi_r(r, theta, phi)     = a(r) Y_lm (theta, phi) exp(-i*2*pi*f*t)
          xi_theta(r, theta, phi) = b(r) dY_lm/dtheta (theta, phi) exp(-i*2*pi*f*t)
          xi_phi(r, theta, phi)   = b(r)/sin(theta) dY_lm/dphi (theta, phi) exp(-i*2*pi*f*t)

        where:

          b(r) = a(r) GM/(R^3*f^2)
        """
        # TODO: we do want to displace the coords_for_observations, but the x,y,z,r below are from the ALSO displaced coords_for_computations
        # if not self._teffext:
            # return coords_for_observations

        x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1))
        theta = np.arccos(z/r)
        phi = np.arctan2(y, x)

        xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
        xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
        xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)

        new_coords = np.zeros(coords_for_observations.shape)
        new_coords[:,0] = coords_for_observations[:,0] + xi_r * np.sin(theta) * np.cos(phi)
        new_coords[:,1] = coords_for_observations[:,1] + xi_r * np.sin(theta) * np.sin(phi)
        new_coords[:,2] = coords_for_observations[:,2] + xi_r * np.cos(theta)

        return new_coords 

def dYdphi(self, m, l, theta, phi):
        return 1j*m*Y(m, l, theta, phi) 

def dYdtheta(self, m, l, theta, phi):
        if abs(m) > l:
            return 0

        # TODO: just a quick hack
        if abs(m+1) > l:
            last_term = 0.0
        else:
            last_term = Y(m+1, l, theta, phi)

        return m/np.tan(theta)*Y(m, l, theta, phi) + np.sqrt((l-m)*(l+m+1))*np.exp(-1j*phi)*last_term 

def sph_harm_lm(l, m, n):
    """ Wrapper around scipy.special.sph_harm. Return spherical harmonic of degree l and order m. """
    phi, theta = util.sph_sample(n)
    phi, theta = np.meshgrid(phi, theta)
    f = sph_harm(m, l, theta, phi)

    return f 

def real_sph_harm(mm, ll, phi, theta):
    """
    The real-valued spherical harmonics.
    """
    if mm > 0:
        ans = (1.0 / np.sqrt(2)) * (ss.sph_harm(mm, ll, phi, theta) + ((-1) ** mm) * ss.sph_harm(-mm, ll, phi, theta))
    elif mm == 0:
        ans = ss.sph_harm(0, ll, phi, theta)
    elif mm < 0:
        ans = (1.0 / (np.sqrt(2) * complex(0.0, 1))) * (
            ss.sph_harm(-mm, ll, phi, theta) - ((-1) ** mm) * ss.sph_harm(mm, ll, phi, theta)
        )

    return ans.real 

def __call__(self, coords):
        from scipy.special import sph_harm

        theta, phi = cart_to_polar_angles(coords)
        if self.m == 0:
            return (sph_harm(self.m, self.l, phi, theta)).real

        vplus = sph_harm(abs(self.m), self.l, phi, theta)
        vminus = sph_harm(-abs(self.m), self.l, phi, theta)
        value = np.sqrt(1/2.0) * (self._posneg*vplus + vminus)

        if self.m < 0:
            return -value.imag
        else:
            return value.real 

def _naive_csh_seismology(l, m, theta, phi):
    """
    Compute the spherical harmonics according to the seismology convention, in a naive way.
    This appears to be equal to the sph_harm function in scipy.special.
    """
    return (lpmv(m, l, np.cos(theta)) * np.exp(1j * m * phi) *
            np.sqrt(((2 * l + 1) * factorial(l - m))
                    /
                    (4 * np.pi * factorial(l + m)))) 

def test_spherharm(self):
        def spherharm(l, m, theta, phi):
            if m > l:
                return np.nan
            return sc.sph_harm(m, l, phi, theta)
        assert_mpmath_equal(spherharm,
                            mpmath.spherharm,
                            [IntArg(0, 100), IntArg(0, 100),
                             Arg(a=0, b=pi), Arg(a=0, b=2*pi)],
                            atol=1e-8, n=6000,
                            dps=150) 

def test_first_harmonics():
    # Test against explicit representations of the first four
    # spherical harmonics which use `theta` as the azimuthal angle,
    # `phi` as the polar angle, and include the Condon-Shortley
    # phase.

    # Notation is Ymn
    def Y00(theta, phi):
        return 0.5*np.sqrt(1/np.pi)

    def Yn11(theta, phi):
        return 0.5*np.sqrt(3/(2*np.pi))*np.exp(-1j*theta)*np.sin(phi)

    def Y01(theta, phi):
        return 0.5*np.sqrt(3/np.pi)*np.cos(phi)

    def Y11(theta, phi):
        return -0.5*np.sqrt(3/(2*np.pi))*np.exp(1j*theta)*np.sin(phi)

    harms = [Y00, Yn11, Y01, Y11]
    m = [0, -1, 0, 1]
    n = [0, 1, 1, 1]

    theta = np.linspace(0, 2*np.pi)
    phi = np.linspace(0, np.pi)
    theta, phi = np.meshgrid(theta, phi)

    for harm, m, n in zip(harms, m, n):
        assert_allclose(sc.sph_harm(m, n, theta, phi),
                        harm(theta, phi),
                        rtol=1e-15, atol=1e-15,
                        err_msg="Y^{}_{} incorrect".format(m, n)) 

def test_sph_harm_ufunc_loop_selection():
    # see https://github.com/scipy/scipy/issues/4895
    dt = np.dtype(np.complex128)
    assert_equal(special.sph_harm(0, 0, 0, 0).dtype, dt)
    assert_equal(special.sph_harm([0], 0, 0, 0).dtype, dt)
    assert_equal(special.sph_harm(0, [0], 0, 0).dtype, dt)
    assert_equal(special.sph_harm(0, 0, [0], 0).dtype, dt)
    assert_equal(special.sph_harm(0, 0, 0, [0]).dtype, dt)
    assert_equal(special.sph_harm([0], [0], [0], [0]).dtype, dt) 

def test_sph_harm():
    # Tests derived from tables in
    # http://en.wikipedia.org/wiki/Table_of_spherical_harmonics
    sh = special.sph_harm
    pi = np.pi
    exp = np.exp
    sqrt = np.sqrt
    sin = np.sin
    cos = np.cos
    assert_array_almost_equal(sh(0,0,0,0),
           0.5/sqrt(pi))
    assert_array_almost_equal(sh(-2,2,0.,pi/4),
           0.25*sqrt(15./(2.*pi)) *
           (sin(pi/4))**2.)
    assert_array_almost_equal(sh(-2,2,0.,pi/2),
           0.25*sqrt(15./(2.*pi)))
    assert_array_almost_equal(sh(2,2,pi,pi/2),
           0.25*sqrt(15/(2.*pi)) *
           exp(0+2.*pi*1j)*sin(pi/2.)**2.)
    assert_array_almost_equal(sh(2,4,pi/4.,pi/3.),
           (3./8.)*sqrt(5./(2.*pi)) *
           exp(0+2.*pi/4.*1j) *
           sin(pi/3.)**2. *
           (7.*cos(pi/3.)**2.-1))
    assert_array_almost_equal(sh(4,4,pi/8.,pi/6.),
           (3./16.)*sqrt(35./(2.*pi)) *
           exp(0+4.*pi/8.*1j)*sin(pi/6.)**4.) 

def test_spherharm(self):
        def spherharm(l, m, theta, phi):
            if m > l:
                return np.nan
            return sc.sph_harm(m, l, phi, theta)
        assert_mpmath_equal(spherharm,
                            mpmath.spherharm,
                            [IntArg(0, 100), IntArg(0, 100),
                             Arg(a=0, b=pi), Arg(a=0, b=2*pi)],
                            atol=1e-8, n=6000,
                            dps=150) 

def an_spherical(q, xq, E_1, E_2, E_0, R, N):
    """
    It computes the analytical solution of the potential of a sphere with
    Nq charges inside.
    Took from Kirkwood (1934).

    Arguments
    ----------
    q  : array, charges.
    xq : array, positions of the charges.
    E_1: float, dielectric constant inside the sphere.
    E_2: float, dielectric constant outside the sphere.
    E_0: float, dielectric constant of vacuum.
    R  : float, radius of the sphere.
    N  : int, number of terms desired in the spherical harmonic expansion.

    Returns
    --------
    PHI: array, reaction potential.
    """

    PHI = numpy.zeros(len(q))
    for K in range(len(q)):
        rho = numpy.sqrt(numpy.sum(xq[K]**2))
        zenit = numpy.arccos(xq[K, 2] / rho)
        azim = numpy.arctan2(xq[K, 1], xq[K, 0])

        phi = 0. + 0. * 1j
        for n in range(N):
            for m in range(-n, n + 1):
                sph1 = special.sph_harm(m, n, zenit, azim)
                cons1 = rho**n / (E_1 * E_0 * R**(2 * n + 1)) * (E_1 - E_2) * (
                    n + 1) / (E_1 * n + E_2 * (n + 1))
                cons2 = 4 * pi / (2 * n + 1)

                for k in range(len(q)):
                    rho_k = numpy.sqrt(numpy.sum(xq[k]**2))
                    zenit_k = numpy.arccos(xq[k, 2] / rho_k)
                    azim_k = numpy.arctan2(xq[k, 1], xq[k, 0])
                    sph2 = numpy.conj(special.sph_harm(m, n, zenit_k, azim_k))
                    phi += cons1 * cons2 * q[K] * rho_k**n * sph1 * sph2

        PHI[K] = numpy.real(phi) / (4 * pi)

    return PHI 

def sph_harm_large(m, n, az, el):
    """Compute spherical harmonics for large orders > 84

    Parameters
    ----------
    m : (int)
        Order of the spherical harmonic. abs(m) <= n

    n : (int)
        Degree of the harmonic, sometimes called l. n >= 0

    az: (float)
        Azimuthal (longitudinal) coordinate [0, 2pi], also called Theta.

    el : (float)
        Elevation (colatitudinal) coordinate [0, pi], also called Phi.

    Returns
    -------
    y_mn : (complex float)
        Complex spherical harmonic of order m and degree n,
        sampled at theta = az, phi = el

    Notes
    -----
    Y_n,m (theta, phi) = ((n - m)! * (2l + 1)) / (4pi * (l + m))^0.5 * exp(i m phi) * P_n^m(cos(theta))
    as per http://dlmf.nist.gov/14.30
    Pmn(z) is the associated Legendre function of the first kind, like scipy.special.lpmv
    scipy.special.lpmn calculates P(0...m 0...n) and its derivative but won't return +inf at high orders
    """
    if _np.all(_np.abs(m) < 84):
        return scy.sph_harm(m, n, az, el)
    else:
        # TODO: confirm that this uses the correct SH definition

        mAbs = _np.abs(m)
        if isinstance(el, _np.ndarray):
            P = _np.empty(el.size)
            for k in range(0, el.size):
                P[k] = scy.lpmn(mAbs, n, _np.cos(el[k]))[0][-1][-1]
        else:
            P = scy.lpmn(mAbs, n, _np.cos(el))[0][-1][-1]
        preFactor1 = _np.sqrt((2 * n + 1) / (4 * _np.pi))
        try:
            preFactor2 = _np.sqrt(fact(n - mAbs) / fact(n + mAbs))
        except OverflowError:  # integer division for very large orders
            preFactor2 = _np.sqrt(fact(n - mAbs) // fact(n + mAbs))

        Y = preFactor1 * preFactor2 * _np.exp(1j * m * az) * P
        if m < 0:
            return _np.conj(Y)
        else:
            return Y 

def plot_sphere_func2(f, grid='Clenshaw-Curtis', beta=None, alpha=None, colormap='jet', fignum=0,  normalize=True):
    # TODO: update this  function now that we have changed the order of axes in f
    import matplotlib.pyplot as plt
    from matplotlib import cm, colors
    from mpl_toolkits.mplot3d import Axes3D
    import numpy as np
    from scipy.special import sph_harm

    if normalize:
        f = (f - np.min(f)) / (np.max(f) - np.min(f))

    if grid == 'Driscoll-Healy':
        b = f.shape[0] // 2
    elif grid == 'Clenshaw-Curtis':
        b = (f.shape[0] - 2) // 2
    elif grid == 'SOFT':
        b = f.shape[0] // 2
    elif grid == 'Gauss-Legendre':
        b = (f.shape[0] - 2) // 2

    if beta is None or alpha is None:
        beta, alpha = meshgrid(b=b, grid_type=grid)

    alpha = np.r_[alpha, alpha[0, :][None, :]]
    beta = np.r_[beta, beta[0, :][None, :]]
    f = np.r_[f, f[0, :][None, :]]

    x = np.sin(beta) * np.cos(alpha)
    y = np.sin(beta) * np.sin(alpha)
    z = np.cos(beta)

    # m, l = 2, 3
    # Calculate the spherical harmonic Y(l,m) and normalize to [0,1]
    # fcolors = sph_harm(m, l, beta, alpha).real
    # fmax, fmin = fcolors.max(), fcolors.min()
    # fcolors = (fcolors - fmin) / (fmax - fmin)
    print(x.shape, f.shape)

    if f.ndim == 2:
        f = cm.gray(f)
        print('2')

    # Set the aspect ratio to 1 so our sphere looks spherical
    fig = plt.figure(figsize=plt.figaspect(1.))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=f ) # cm.gray(f))
    # Turn off the axis planes
    ax.set_axis_off()
    plt.show() 

def sht_matrix(N, azi, colat, weights=None):
    r"""Matrix of spherical harmonics up to order N for given angles.

    Computes a matrix of spherical harmonics up to order :math:`N`
    for the given angles/grid.

    .. math::

        \mathbf{Y} = \left[ \begin{array}{cccccc}
        Y_0^0(\theta[0], \phi[0]) & Y_1^{-1}(\theta[0], \phi[0]) & Y_1^0(\theta[0], \phi[0]) & Y_1^1(\theta[0], \phi[0]) & \dots & Y_N^N(\theta[0], \phi[0])  \\
        Y_0^0(\theta[1], \phi[1]) & Y_1^{-1}(\theta[1], \phi[1]) & Y_1^0(\theta[1], \phi[1]) & Y_1^1(\theta[1], \phi[1]) & \dots & Y_N^N(\theta[1], \phi[1])  \\
        \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
        Y_0^0(\theta[Q-1], \phi[Q-1]) & Y_1^{-1}(\theta[Q-1], \phi[Q-1]) & Y_1^0(\theta[Q-1], \phi[Q-1]) & Y_1^1(\theta[Q-1], \phi[Q-1]) & \dots & Y_N^N(\theta[Q-1], \phi[Q-1])
        \end{array} \right]

    where

    .. math::

        Y_n^m(\theta, \phi) = \sqrt{\frac{2n + 1}{4 \pi} \frac{(n-m)!}{(n+m)!}} P_n^m(\cos \theta) e^{i m \phi}


    (Note: :math:`\mathbf{Y}` is interpreted as the inverse transform (or synthesis)
    matrix in examples and documentation.)

    Parameters
    ----------
    N : int
        Maximum order.
    azi : (Q,) array_like
        Azimuth.
    colat : (Q,) array_like
        Colatitude.
    weights : (Q,) array_like, optional
        Quadrature weights.

    Returns
    -------
    Ymn : (Q, (N+1)**2) numpy.ndarray
        Matrix of spherical harmonics.

    """
    azi = util.asarray_1d(azi)
    colat = util.asarray_1d(colat)
    if azi.ndim == 0:
        Q = 1
    else:
        Q = len(azi)
    if weights is None:
        weights = np.ones(Q)
    Ymn = np.zeros([Q, (N+1)**2], dtype=complex)
    i = 0
    for n in range(N+1):
        for m in range(-n, n+1):
            Ymn[:, i] = weights * special.sph_harm(m, n, azi, colat)
            i += 1
    return Ymn