Python scipy.special 模块，j1() 实例源码

def smooth_fft_vsini(dv,spec,sigma):
    # The Fourier coordinate
    ss = rfftfreq(len(spec), d=dv)

    # Make the fourier space taper
    ss[0] = 0.01 #junk so we don't get a divide by zero error
    ub = 2. * np.pi * sigma * ss
    sb = j1(ub) / ub - 3 * np.cos(ub) / (2 * ub ** 2) + 3. * np.sin(ub) / (2 * ub ** 3)
    #set zeroth frequency to 1 separately (DC term)
    sb[0] = 1.

    # Fourier transform the spectrum
    FF = np.fft.rfft(spec)

    # Multiply in fourier space
    FF_tap = FF * sb

    # Fourier transform back
    spec_conv = np.fft.irfft(FF_tap)
    return spec_conv

def jinc(r):
    ''' The Jinc function.

    Args:
        r (`number`): radial distance.

    Returns:
        `float`: the value of j1(x)/x for x != 0, 0.5 at 0.

    '''
    if r == 0:
        return 0.5
    else:
        return j1(r) / r

def LambDipole(model, U=.01,R = 1.):

    """ Generate Lamb's dipole vorticity field.

       Parameters
       -----------
        U: float
                Translation speed (dipole's strength).
        R: float
                Diple's radius.

       Return
       -------
        q: array of floats
              Vorticity (physical space).

    """

    N = model.nx
    x, y = model.x, model.y
    x0,y0 = x[N//2,N//2],y[N//2,N//2]

    r = np.sqrt( (x-x0)**2 + (y-y0)**2 )
    s = np.zeros_like(r)

    for i in range(N):
        for j in range(N):
            if r[i,j] == 0.:
                s[i,j] = 0.
            else:
                s[i,j] = (y[i,j]-y0)/r[i,j]

    lam = (3.8317)/R
    C = -(2.*U*lam)/(special.j0(lam*R))
    q = np.zeros_like(r)
    q[r<=R] = C*special.j1(lam*r[r<=R])*s[r<=R]

    return q

def airy(th, B, lam=None):
    """
    Return the visibility value (unsquared), given
    - th the angular diameter in radian
    - B the baseline in m (alternatively in m/lambda)
    - lam the wavelength (alternatively None if B is already given as m/lambda)
    """
    if lam is None:
        x = np.pi*th*B
    else:
        x = np.pi*th*B/lam
    return 2*airyJ1(x)/x

def lanczos(n, fc=0.02):
    """
    Compute the coefficients of a Lanczos window

    Parameters
    ----------
    n : int
        Number of points in the output window, must be an odd integer.
    fc : float
        Cutoff frequency

    Returns
    -------
    w : ndarray
        The weights associated to the boxcar window
    """
    if not isinstance(n, int):
        try:
            n = int(n)
        except:
            TypeError, "n must be an integer"
    if not n % 2 == 1:
        raise ValueError, "n must an odd integer"
    dim = 1
    if dim ==  1:
        k = np.arange(- n / 2 + 1, n / 2 + 1)
        # k = np.arange(0, n + 1)
        # w = (np.sin(2 * pi * fc * k) / (pi * k) *
        #     np.sin(pi * k / n) / (pi * k / n))
        w = (np.sin(2. * np.pi * fc * k) / (np.pi * k) *
             np.sin(np.pi * k / (n / 2.)) / (np.pi * k / (n / 2.)))
        # Particular case where k=0
        w[n / 2] = 2. * fc
    elif dim == 2:
        #TODO: Test this bidimensional window
        fcx, fcy = fc
        nx, ny = n
        # Grid definition according to the number of weights
        kx, ky = np.meshgrid(np.arange(-nx, nx + 1), np.arange(-ny, ny + 1), indexing='ij')
        # Computation of the response weight on the grid
        z = np.sqrt((fcx * kx) ** 2 + (fcy * ky) ** 2)
        w_rect = 1 / z * fcx * fcy * spec.j1(2 * np.pi * z)
        w = (w_rect * np.sin(np.pi * kx / nx) / (np.pi * kx / nx) *
             np.sin(np.pi * ky / ny) / (np.pi * ky / ny))
        # Particular case where z=0
        w[nx, :] = (w_rect[nx, :] * 1. / (np.pi * ky[nx, :] / ny) *
                     np.sin(np.pi * ky[nx, :] / ny))
        w[:, ny] = (w_rect[:, ny] * 1. / (np.pi * kx[:, ny] / nx) *
                     np.sin(np.pi * kx[:, ny] / nx))
        w[nx, ny] = np.pi * fcx * fcy
    return w


# Define a list of available windows

def quad(sPJ0r, sPJ0i, sPJ1r, sPJ1i, sPJ0br, sPJ0bi, ab, off, factAng, iinp):
    """Quadrature for Hankel transform.

    This is the kernel of the QUAD method, used for the Hankel transforms
    `hquad` and `hqwe` (where the integral is not suited for QWE).

    """

    # Define the quadrature kernels
    def quad0(klambd, sPJ, sPJb, koff, kang):
        """Quadrature for J0."""
        tP0 = sPJ(np.log10(klambd)) + kang*sPJb(np.log10(klambd))
        return tP0*special.j0(koff*klambd)

    def quad1(klambd, sPJ, ab, koff, kang):
        """Quadrature for J1."""
        tP1 = kang*sPJ(np.log10(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)

    # Carry out quadrature of J0
    iargs = (sPJ0r, sPJ0br, off, factAng)
    fr0 = integrate.quad(quad0, args=iargs, full_output=1, **iinp)
    iargs = (sPJ0i, sPJ0bi, off, factAng)
    fi0 = integrate.quad(quad0, args=iargs, full_output=1, **iinp)

    # Carry out quadrature of J1
    iargs = (sPJ1r, ab, off, factAng)
    fr1 = integrate.quad(quad1, args=iargs, full_output=1, **iinp)
    iargs = (sPJ1i, ab, off, factAng)
    fi1 = integrate.quad(quad1, args=iargs, full_output=1, **iinp)

    # If there is a fourth output from QUAD, it means it did not converge
    if np.any(np.array([len(fr0), len(fi0), len(fr1), len(fi1)]) > 3):
        conv = False
    else:
        conv = True

    # Collect the results
    return fr0[0] + fr1[0] + 1j*(fi0[0] + fi1[0]), conv