transform.py 文件源码

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