Python special.j1方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def _encircled_energy_core(mtf_data, radius, nu_p, dx, dy):
    """Core computation of encircled energy, based on Baliga 1988.

    Parameters
    ----------
    mtf_data : `numpy.ndarray`
        unaliased MTF data
    radius : `float`
        radius of "detector"
    nu_p : `numpy.ndarray`
        radial spatial frequencies
    dx : `float`
        x frequency delta
    dy : `float`
        y frequency delta

    Returns
    -------
    `float`
        encircled energy for given radius

    """
    integration_fourier = special.j1(2 * e.pi * radius * nu_p) / nu_p
    dat = mtf_data * integration_fourier
    return radius * dat.sum() * dx * dy 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def _analytical_encircled_energy(fno, wavelength, points):
    """Compute the analytical encircled energy for a diffraction limited circular aperture.

    Parameters
    ----------
    fno : `float`
        F/#
    wavelength : `float`
        wavelength of light
    points : `numpy.ndarray`
        radii of "detector"

    Returns
    -------
    `numpy.ndarray`
        encircled energy values

    """
    p = points * e.pi / fno / wavelength
    return 1 - special.j0(p)**2 - special.j1(p)**2 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def spectral_density(self, k):  # noqa: D102
        k = np.array(k, dtype=np.double)
        res = np.empty_like(k)
        kl = k * self.len_scale
        kl_gz = kl > 0
        # for k=0 we calculate the limit by hand
        if self.dim == 1:
            res[kl_gz] = (1.0 - np.cos(kl[kl_gz])) / (
                np.pi * k[kl_gz] * kl[kl_gz]
            )
            res[np.logical_not(kl_gz)] = self.len_scale / 2.0 / np.pi
        elif self.dim == 2:
            res[kl_gz] = sps.j1(kl[kl_gz] / 2.0) ** 2 / np.pi / k[kl_gz] ** 2
            res[np.logical_not(kl_gz)] = self.len_scale ** 2 / 16.0 / np.pi
        else:
            res[kl_gz] = -(
                12 * kl[kl_gz] * np.sin(kl[kl_gz])
                + (12 - 3 * kl[kl_gz] ** 2) * np.cos(kl[kl_gz])
                - 3 * kl[kl_gz] ** 2
                - 12
            ) / (2 * np.pi ** 2 * kl[kl_gz] ** 3 * k[kl_gz] ** 3)
            res[np.logical_not(kl_gz)] = (
                self.len_scale ** 3 / 48.0 / np.pi ** 2
            )
        return res 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def airy_disk(x):   
    return np.power( 2.*j1(x)/(x) , 2)


# generate a (2*size+1,2*size+1) convolution kernel with "radii" scale
# where the function above is assumed to have "radius" one
# scale is a 3-vector for RGB 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def test_j1(self):
        assert_equal(cephes.j1(0),0.0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def test_j1(self):
        # See comment in test_j0
        assert_mpmath_equal(sc.j1,
                            mpmath.j1,
                            [Arg(-1e3, 1e3)])
        assert_mpmath_equal(sc.j1,
                            mpmath.j1,
                            [Arg(-1e8, 1e8)],
                            rtol=1e-5) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def perpendicular_attenuation(self, q, tau, diameter):
        "Implements the finite time Callaghan model for cylinders"
        radius = diameter / 2.
        alpha = self.alpha
        q_argument = 2 * np.pi * q * radius
        q_argument_2 = q_argument ** 2
        res = np.zeros_like(q)

        J = special.j1(q_argument) ** 2
        for k in range(0, self.alpha.shape[0]):
            alpha2 = alpha[k, 0] ** 2
            update = (
                4 * np.exp(-alpha2 * self.diffusion_perpendicular *
                           tau / radius ** 2) *
                q_argument_2 /
                (q_argument_2 - alpha2) ** 2 * J
            )
            res += update

        for m in range(1, self.alpha.shape[1]):
            J = special.jvp(m, q_argument, 1)
            q_argument_J = (q_argument * J) ** 2
            for k in range(self.alpha.shape[0]):
                alpha2 = self.alpha[k, m] ** 2
                update = (
                    8 * np.exp(-alpha2 * self.diffusion_perpendicular *
                               tau / radius ** 2) *
                    alpha2 / (alpha2 - m ** 2) *
                    q_argument_J /
                    (q_argument_2 - alpha2) ** 2
                )
                res += update
        return res 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def test_radial_reduction():
    """radial reduction example"""
    x = np.linspace(-2,2, 64)
    y = x[:, None]
    x = x[None, :]
    R = np.sqrt(x**2+y**2)

    def airy(r, sigma):
        from scipy.special import j1
        r = r / sigma * np.sqrt(2)
        a = (2*j1(r)/r)**2
        a[r==0] = 1
        return a
    def gauss(r, sigma):
        return np.exp(-(r/sigma)**2)

    distribution = np.random.choice([gauss, airy])(R, 0.3)
    sample = np.random.poisson(distribution*200+10).astype(np.float)

    #is this an airy or gaussian function? hard to tell with all this noise!
    plt.imshow(sample, interpolation='nearest', cmap='gray')
    plt.show()
    #radial reduction to the rescue!
    #if we are sampling an airy function, you will see a small but significant rise around x=1
    g = group_by(np.round(R, 5).flatten())
    plt.errorbar(
        g.unique,
        g.mean(sample.flatten())[1],
        g.std (sample.flatten())[1] / np.sqrt(g.count))
    plt.xlim(0, 2)
    plt.show() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def sears_lift_sin_gust(w0, L, Uinf, chord, tv):
    """
    Returns the lift coefficient for a sinusoidal gust (see set_gust.sin) as
    the imaginary part of the CL complex function defined below. The input gust
    must be the imaginary part of

    .. math::    wgust = w0*\exp(1.0j*C*(Ux*S.time[tt] - xcoord) )

    with:

    .. math:: C=2\pi/L

    and ``xcoord=0`` at the aerofoil half-chord.
    """

    # reduced frequency
    kg = np.pi * chord / L
    # Theo's funciton
    Ctheo = theo_fun(kg)
    # Sear's function
    J0, J1 = scsp.j0(kg), scsp.j1(kg)
    S = (J0 - 1.0j * J1) * Ctheo + 1.0j * J1

    phase = np.angle(S)
    CL = 2. * np.pi * w0 / Uinf * np.abs(S) * np.sin(2. * np.pi * Uinf / L * tv + phase)

    return CL 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import j1 [as 别名]
def quad(sPJ0r, sPJ0i, sPJ1r, sPJ1i, sPJ0br, sPJ0bi, ab, off, ang_fact, iinp):
    r"""Quadrature for Hankel transform.

    This is the kernel of the QUAD method, used for the Hankel transforms
    :func:`hankel_quad` and :func:`hankel_qwe` (where the integral is not
    suited for QWE).

    """

    # Define the quadrature kernels
    def quad_PJ0(klambd, sPJ0, koff):
        r"""Quadrature for PJ0."""
        return sPJ0(np.log(klambd))*special.j0(koff*klambd)

    def quad_PJ1(klambd, sPJ1, ab, koff, kang):
        r"""Quadrature for PJ1."""

        tP1 = kang*sPJ1(np.log(klambd))
        if ab in [11, 12, 21, 22, 14, 24, 15, 25]:  # Because of J2
            # J2(kr) = 2/(kr)*J1(kr) - J0(kr)
            tP1 /= koff

        return tP1*special.j1(koff*klambd)

    def quad_PJ0b(klambd, sPJ0b, koff, kang):
        r"""Quadrature for PJ0b."""
        return kang*sPJ0b(np.log(klambd))*special.j0(koff*klambd)

    # Pre-allocate output
    conv = True
    out = np.array(0.0+0.0j)

    # Carry out quadrature for required kernels
    iinp['full_output'] = 1

    if sPJ0r is not None:
        re = integrate.quad(quad_PJ0, args=(sPJ0r, off), **iinp)
        im = integrate.quad(quad_PJ0, args=(sPJ0i, off), **iinp)
        out += re[0] + 1j*im[0]
        # If there is a fourth output from QUAD, it means it did not converge
        if (len(re) or len(im)) > 3:
            conv = False

    if sPJ1r is not None:
        re = integrate.quad(quad_PJ1, args=(sPJ1r, ab, off, ang_fact), **iinp)
        im = integrate.quad(quad_PJ1, args=(sPJ1i, ab, off, ang_fact), **iinp)
        out += re[0] + 1j*im[0]
        # If there is a fourth output from QUAD, it means it did not converge
        if (len(re) or len(im)) > 3:
            conv = False

    if sPJ0br is not None:
        re = integrate.quad(quad_PJ0b, args=(sPJ0br, off, ang_fact), **iinp)
        im = integrate.quad(quad_PJ0b, args=(sPJ0bi, off, ang_fact), **iinp)
        out += re[0] + 1j*im[0]
        # If there is a fourth output from QUAD, it means it did not converge
        if (len(re) or len(im)) > 3:
            conv = False

    # Collect the results
    return out, conv