Python scipy.special.j0() Examples

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 

def test_j0(self):
        assert_equal(cephes.j0(0),1.0) 

def test_j0(self):
        oz = special.j0(.1)
        ozr = special.jn(0,.1)
        assert_almost_equal(oz,ozr,8) 

def test_j0(self):
        # The Bessel function at large arguments is j0(x) ~ cos(x + phi)/sqrt(x)
        # and at large arguments the phase of the cosine loses precision.
        #
        # This is numerically expected behavior, so we compare only up to
        # 1e8 = 1e15 * 1e-7
        assert_mpmath_equal(sc.j0,
                            mpmath.j0,
                            [Arg(-1e3, 1e3)])
        assert_mpmath_equal(sc.j0,
                            mpmath.j0,
                            [Arg(-1e8, 1e8)],
                            rtol=1e-5) 

def test_j0(self):
        assert_equal(cephes.j0(0),1.0) 

def test_j0(self):
        oz = special.j0(.1)
        ozr = special.jn(0,.1)
        assert_almost_equal(oz,ozr,8) 

def test_j0(self):
        # The Bessel function at large arguments is j0(x) ~ cos(x + phi)/sqrt(x)
        # and at large arguments the phase of the cosine loses precision.
        #
        # This is numerically expected behavior, so we compare only up to
        # 1e8 = 1e15 * 1e-7
        assert_mpmath_equal(sc.j0,
                            mpmath.j0,
                            [Arg(-1e3, 1e3)])
        assert_mpmath_equal(sc.j0,
                            mpmath.j0,
                            [Arg(-1e8, 1e8)],
                            rtol=1e-5) 

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 

def _hankel2_0(x):
    """Wrapper for Hankel function of the second type using fast versions
       of the Bessel functions of first/second kind in scipy"""
    return _special.j0(x) - 1j * _special.y0(x) 

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