specialpython_Python special.j1方法代码示例

# 需要导入模块: 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