Python special.sph_harm方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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.) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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.) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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)))) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
def dYdphi(self, m, l, theta, phi):
        return 1j*m*Y(m, l, theta, phi) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import sph_harm [as 别名]
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