basic.py

#
# Author:  Travis Oliphant, 2002
#

from numpy import pi, asarray, floor, isscalar, iscomplex, real, imag, sqrt, \
        where, mgrid, cos, sin, exp, place, seterr, issubdtype, extract, \
        complexfloating, less, vectorize, inexact, nan, zeros, sometrue
from _cephes import ellipkm1, mathieu_a, mathieu_b, iv, jv, gamma, psi, zeta, \
        hankel1, hankel2, yv, kv, gammaln, errprint, ndtri
import types
import specfun
import orthogonal

__all__ = ['agm', 'ai_zeros', 'assoc_laguerre', 'bei_zeros', 'beip_zeros',
           'ber_zeros', 'bernoulli', 'berp_zeros', 'bessel_diff_formula',
           'bi_zeros', 'digamma', 'diric', 'ellipk', 'erf_zeros', 'erfcinv',
           'erfinv', 'errprint', 'euler', 'fresnel_zeros',
           'fresnelc_zeros', 'fresnels_zeros', 'gamma', 'gammaln', 'h1vp',
           'h2vp', 'hankel1', 'hankel2', 'hyp0f1', 'iv', 'ivp', 'jn_zeros',
           'jnjnp_zeros', 'jnp_zeros', 'jnyn_zeros', 'jv', 'jvp', 'kei_zeros',
           'keip_zeros', 'kelvin_zeros', 'ker_zeros', 'kerp_zeros', 'kv',
           'kvp', 'lmbda', 'lpmn', 'lpn', 'lqmn', 'lqn', 'mathieu_a',
           'mathieu_b', 'mathieu_even_coef', 'mathieu_odd_coef', 'ndtri',
           'obl_cv_seq', 'pbdn_seq', 'pbdv_seq', 'pbvv_seq',
           'polygamma', 'pro_cv_seq', 'psi', 'riccati_jn', 'riccati_yn',
           'sinc', 'sph_harm', 'sph_in', 'sph_inkn',
           'sph_jn', 'sph_jnyn', 'sph_kn', 'sph_yn', 'y0_zeros', 'y1_zeros',
           'y1p_zeros', 'yn_zeros', 'ynp_zeros', 'yv', 'yvp', 'zeta']


def sinc(x):
    """Returns sin(pi*x)/(pi*x) at all points of array x.
    """
    w = pi * asarray(x)
    # w might contain 0, and so temporarily turn off warnings
    # while calculating sin(w)/w.
    old_settings = seterr(all='ignore')
    s = sin(w) / w
    seterr(**old_settings)
    return where(x==0, 1.0, s)

def diric(x,n):
    """Returns the periodic sinc function also called the dirichlet function:

    diric(x) = sin(x *n / 2) / (n sin(x / 2))

    where n is a positive integer.
    """
    x,n = asarray(x), asarray(n)
    n = asarray(n + (x-x))
    x = asarray(x + (n-n))
    if issubdtype(x.dtype, inexact):
        ytype = x.dtype
    else:
        ytype = float
    y = zeros(x.shape,ytype)

    mask1 = (n <= 0) | (n <> floor(n))
    place(y,mask1,nan)

    z = asarray(x / 2.0 / pi)
    mask2 = (1-mask1) & (z == floor(z))
    zsub = extract(mask2,z)
    nsub = extract(mask2,n)
    place(y,mask2,pow(-1,zsub*(nsub-1)))

    mask = (1-mask1) & (1-mask2)
    xsub = extract(mask,x)
    nsub = extract(mask,n)
    place(y,mask,sin(nsub*xsub/2.0)/(nsub*sin(xsub/2.0)))
    return y



def jnjnp_zeros(nt):
    """Compute nt (<=1200) zeros of the bessel functions Jn and Jn'
    and arange them in order of their magnitudes.

    Returns
    -------
    zo[l-1] : ndarray
        Value of the lth zero of of Jn(x) and Jn'(x). Of length `nt`.
    n[l-1] : ndarray
        Order of the Jn(x) or Jn'(x) associated with lth zero. Of length `nt`.
    m[l-1] : ndarray
        Serial number of the zeros of Jn(x) or Jn'(x) associated
        with lth zero. Of length `nt`.
    t[l-1] : ndarray
        0 if lth zero in zo is zero of Jn(x), 1 if it is a zero of Jn'(x). Of
        length `nt`.

    See Also
    --------
    jn_zeros, jnp_zeros : to get separated arrays of zeros.
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt>1200):
        raise ValueError("Number must be integer <= 1200.")
    nt = int(nt)
    n,m,t,zo = specfun.jdzo(nt)
    return zo[1:nt+1],n[:nt],m[:nt],t[:nt]

def jnyn_zeros(n,nt):
    """Compute nt zeros of the Bessel functions Jn(x), Jn'(x), Yn(x), and
    Yn'(x), respectively. Returns 4 arrays of length nt.

    See jn_zeros, jnp_zeros, yn_zeros, ynp_zeros to get separate arrays.
    """
    if not (isscalar(nt) and isscalar(n)):
        raise ValueError("Arguments must be scalars.")
    if (floor(n)!=n) or (floor(nt)!=nt):
        raise ValueError("Arguments must be integers.")
    if (nt <=0):
        raise ValueError("nt > 0")
    return specfun.jyzo(abs(n),nt)

def jn_zeros(n,nt):
    """Compute nt zeros of the Bessel function Jn(x).
    """
    return jnyn_zeros(n,nt)[0]
def jnp_zeros(n,nt):
    """Compute nt zeros of the Bessel function Jn'(x).
    """
    return jnyn_zeros(n,nt)[1]
def yn_zeros(n,nt):
    """Compute nt zeros of the Bessel function Yn(x).
    """
    return jnyn_zeros(n,nt)[2]
def ynp_zeros(n,nt):
    """Compute nt zeros of the Bessel function Yn'(x).
    """
    return jnyn_zeros(n,nt)[3]

def y0_zeros(nt,complex=0):
    """Returns nt (complex or real) zeros of Y0(z), z0, and the value
    of Y0'(z0) = -Y1(z0) at each zero.
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt <=0):
        raise ValueError("Arguments must be scalar positive integer.")
    kf = 0
    kc = (complex != 1)
    return specfun.cyzo(nt,kf,kc)

def y1_zeros(nt,complex=0):
    """Returns nt (complex or real) zeros of Y1(z), z1, and the value
    of Y1'(z1) = Y0(z1) at each zero.
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt <=0):
        raise ValueError("Arguments must be scalar positive integer.")
    kf = 1
    kc = (complex != 1)
    return specfun.cyzo(nt,kf,kc)

def y1p_zeros(nt,complex=0):
    """Returns nt (complex or real) zeros of Y1'(z), z1', and the value
    of Y1(z1') at each zero.
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt <=0):
        raise ValueError("Arguments must be scalar positive integer.")
    kf = 2
    kc = (complex != 1)
    return specfun.cyzo(nt,kf,kc)

def bessel_diff_formula(v, z, n, L, phase):
    # from AMS55.
    # L(v,z) = J(v,z), Y(v,z), H1(v,z), H2(v,z), phase = -1
    # L(v,z) = I(v,z) or exp(v*pi*i)K(v,z), phase = 1
    # For K, you can pull out the exp((v-k)*pi*i) into the caller
    p = 1.0
    s = L(v-n, z)
    for i in xrange(1, n+1):
        p = phase * (p * (n-i+1)) / i   # = choose(k, i)
        s += p*L(v-n + i*2, z)
    return s / (2.**n)

def jvp(v,z,n=1):
    """Return the nth derivative of Jv(z) with respect to z.
    """
    if not isinstance(n,types.IntType) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if n == 0:
        return jv(v,z)
    else:
        return bessel_diff_formula(v, z, n, jv, -1)
#        return (jvp(v-1,z,n-1) - jvp(v+1,z,n-1))/2.0

def yvp(v,z,n=1):
    """Return the nth derivative of Yv(z) with respect to z.
    """
    if not isinstance(n,types.IntType) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if n == 0:
        return yv(v,z)
    else:
        return bessel_diff_formula(v, z, n, yv, -1)
#        return (yvp(v-1,z,n-1) - yvp(v+1,z,n-1))/2.0

def kvp(v,z,n=1):
    """Return the nth derivative of Kv(z) with respect to z.
    """
    if not isinstance(n,types.IntType) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if n == 0:
        return kv(v,z)
    else:
        return (-1)**n * bessel_diff_formula(v, z, n, kv, 1)

def ivp(v,z,n=1):
    """Return the nth derivative of Iv(z) with respect to z.
    """
    if not isinstance(n,types.IntType) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if n == 0:
        return iv(v,z)
    else:
        return bessel_diff_formula(v, z, n, iv, 1)

def h1vp(v,z,n=1):
    """Return the nth derivative of H1v(z) with respect to z.
    """
    if not isinstance(n,types.IntType) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if n == 0:
        return hankel1(v,z)
    else:
        return bessel_diff_formula(v, z, n, hankel1, -1)
#        return (h1vp(v-1,z,n-1) - h1vp(v+1,z,n-1))/2.0

def h2vp(v,z,n=1):
    """Return the nth derivative of H2v(z) with respect to z.
    """
    if not isinstance(n,types.IntType) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if n == 0:
        return hankel2(v,z)
    else:
        return bessel_diff_formula(v, z, n, hankel2, -1)
#        return (h2vp(v-1,z,n-1) - h2vp(v+1,z,n-1))/2.0

def sph_jn(n,z):
    """Compute the spherical Bessel function jn(z) and its derivative for
    all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n < 1): n1 = 1
    else: n1 = n
    if iscomplex(z):
        nm,jn,jnp,yn,ynp = specfun.csphjy(n1,z)
    else:
        nm,jn,jnp = specfun.sphj(n1,z)
    return jn[:(n+1)], jnp[:(n+1)]

def sph_yn(n,z):
    """Compute the spherical Bessel function yn(z) and its derivative for
    all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n < 1): n1 = 1
    else: n1 = n
    if iscomplex(z) or less(z,0):
        nm,jn,jnp,yn,ynp = specfun.csphjy(n1,z)
    else:
        nm,yn,ynp = specfun.sphy(n1,z)
    return yn[:(n+1)], ynp[:(n+1)]

def sph_jnyn(n,z):
    """Compute the spherical Bessel functions, jn(z) and yn(z) and their
    derivatives for all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n < 1): n1 = 1
    else: n1 = n
    if iscomplex(z) or less(z,0):
        nm,jn,jnp,yn,ynp = specfun.csphjy(n1,z)
    else:
        nm,yn,ynp = specfun.sphy(n1,z)
        nm,jn,jnp = specfun.sphj(n1,z)
    return jn[:(n+1)],jnp[:(n+1)],yn[:(n+1)],ynp[:(n+1)]

def sph_in(n,z):
    """Compute the spherical Bessel function in(z) and its derivative for
    all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n < 1): n1 = 1
    else: n1 = n
    if iscomplex(z):
        nm,In,Inp,kn,knp = specfun.csphik(n1,z)
    else:
        nm,In,Inp = specfun.sphi(n1,z)
    return In[:(n+1)], Inp[:(n+1)]

def sph_kn(n,z):
    """Compute the spherical Bessel function kn(z) and its derivative for
    all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n < 1): n1 = 1
    else: n1 = n
    if iscomplex(z) or less(z,0):
        nm,In,Inp,kn,knp = specfun.csphik(n1,z)
    else:
        nm,kn,knp = specfun.sphk(n1,z)
    return kn[:(n+1)], knp[:(n+1)]

def sph_inkn(n,z):
    """Compute the spherical Bessel functions, in(z) and kn(z) and their
    derivatives for all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if iscomplex(z) or less(z,0):
        nm,In,Inp,kn,knp = specfun.csphik(n,z)
    else:
        nm,In,Inp = specfun.sphi(n,z)
        nm,kn,knp = specfun.sphk(n,z)
    return In,Inp,kn,knp

def riccati_jn(n,x):
    """Compute the Ricatti-Bessel function of the first kind and its
    derivative for all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(x)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n == 0): n1 = 1
    else: n1 = n
    nm,jn,jnp = specfun.rctj(n1,x)
    return jn[:(n+1)],jnp[:(n+1)]

def riccati_yn(n,x):
    """Compute the Ricatti-Bessel function of the second kind and its
    derivative for all orders up to and including n.
    """
    if not (isscalar(n) and isscalar(x)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n == 0): n1 = 1
    else: n1 = n
    nm,jn,jnp = specfun.rcty(n1,x)
    return jn[:(n+1)],jnp[:(n+1)]

def _sph_harmonic(m,n,theta,phi):
    """Compute spherical harmonics.

    This is a ufunc and may take scalar or array arguments like any
    other ufunc.  The inputs will be broadcasted against each other.

    Parameters
    ----------
    m : int
       |m| <= n; the order of the harmonic.
    n : int
       where `n` >= 0; the degree of the harmonic.  This is often called
       ``l`` (lower case L) in descriptions of spherical harmonics.
    theta : float
       [0, 2*pi]; the azimuthal (longitudinal) coordinate.
    phi : float
       [0, pi]; the polar (colatitudinal) coordinate.

    Returns
    -------
    y_mn : complex float
       The harmonic $Y^m_n$ sampled at `theta` and `phi`

    Notes
    -----
    There are different conventions for the meaning of input arguments
    `theta` and `phi`.  We take `theta` to be the azimuthal angle and
    `phi` to be the polar angle.  It is common to see the opposite
    convention - that is `theta` as the polar angle and `phi` as the
    azimuthal angle.
    """
    x = cos(phi)
    m,n = int(m), int(n)
    Pmn,Pmn_deriv = lpmn(m,n,x)
    # Legendre call generates all orders up to m and degrees up to n
    val = Pmn[-1, -1]
    val *= sqrt((2*n+1)/4.0/pi)
    val *= exp(0.5*(gammaln(n-m+1)-gammaln(n+m+1)))
    val *= exp(1j*m*theta)
    return val

sph_harm = vectorize(_sph_harmonic,'D')

def erfinv(y):
    return ndtri((y+1)/2.0)/sqrt(2)

def erfcinv(y):
    return ndtri((2-y)/2.0)/sqrt(2)

def erf_zeros(nt):
    """Compute nt complex zeros of the error function erf(z).
    """
    if (floor(nt)!=nt) or (nt<=0) or not isscalar(nt):
        raise ValueError("Argument must be positive scalar integer.")
    return specfun.cerzo(nt)

def fresnelc_zeros(nt):
    """Compute nt complex zeros of the cosine fresnel integral C(z).
    """
    if (floor(nt)!=nt) or (nt<=0) or not isscalar(nt):
        raise ValueError("Argument must be positive scalar integer.")
    return specfun.fcszo(1,nt)

def fresnels_zeros(nt):
    """Compute nt complex zeros of the sine fresnel integral S(z).
    """
    if (floor(nt)!=nt) or (nt<=0) or not isscalar(nt):
        raise ValueError("Argument must be positive scalar integer.")
    return specfun.fcszo(2,nt)

def fresnel_zeros(nt):
    """Compute nt complex zeros of the sine and cosine fresnel integrals
    S(z) and C(z).
    """
    if (floor(nt)!=nt) or (nt<=0) or not isscalar(nt):
        raise ValueError("Argument must be positive scalar integer.")
    return specfun.fcszo(2,nt), specfun.fcszo(1,nt)

def hyp0f1(v,z):
    """Confluent hypergeometric limit function 0F1.
    Limit as q->infinity of 1F1(q;a;z/q)
    """
    z = asarray(z)
    if issubdtype(z.dtype, complexfloating):
        arg = 2*sqrt(abs(z))
        num = where(z>=0, iv(v-1,arg), jv(v-1,arg))
        den = abs(z)**((v-1.0)/2)
    else:
        num = iv(v-1,2*sqrt(z))
        den = z**((v-1.0)/2.0)
    num *= gamma(v)
    return where(z==0,1.0,num/ asarray(den))

def assoc_laguerre(x,n,k=0.0):
    return orthogonal.eval_genlaguerre(n, k, x)

digamma = psi

def polygamma(n, x):
    """Polygamma function which is the nth derivative of the digamma (psi)
    function."""
    n, x = asarray(n), asarray(x)
    cond = (n==0)
    fac2 = (-1.0)**(n+1) * gamma(n+1.0) * zeta(n+1,x)
    if sometrue(cond,axis=0):
        return where(cond, psi(x), fac2)
    return fac2

def mathieu_even_coef(m,q):
    """Compute expansion coefficients for even mathieu functions and
    modified mathieu functions.
    """
    if not (isscalar(m) and isscalar(q)):
        raise ValueError("m and q must be scalars.")
    if (q < 0):
        raise ValueError("q >=0")
    if (m != floor(m)) or (m<0):
        raise ValueError("m must be an integer >=0.")

    if (q <= 1):
        qm = 7.5+56.1*sqrt(q)-134.7*q+90.7*sqrt(q)*q
    else:
        qm=17.0+3.1*sqrt(q)-.126*q+.0037*sqrt(q)*q
    km = int(qm+0.5*m)
    if km > 251:
        print "Warning, too many predicted coefficients."
    kd = 1
    m = int(floor(m))
    if m % 2:
        kd = 2

    a = mathieu_a(m,q)
    fc = specfun.fcoef(kd,m,q,a)
    return fc[:km]

def mathieu_odd_coef(m,q):
    """Compute expansion coefficients for even mathieu functions and
    modified mathieu functions.
    """
    if not (isscalar(m) and isscalar(q)):
        raise ValueError("m and q must be scalars.")
    if (q < 0):
        raise ValueError("q >=0")
    if (m != floor(m)) or (m<=0):
        raise ValueError("m must be an integer > 0")

    if (q <= 1):
        qm = 7.5+56.1*sqrt(q)-134.7*q+90.7*sqrt(q)*q
    else:
        qm=17.0+3.1*sqrt(q)-.126*q+.0037*sqrt(q)*q
    km = int(qm+0.5*m)
    if km > 251:
        print "Warning, too many predicted coefficients."
    kd = 4
    m = int(floor(m))
    if m % 2:
        kd = 3

    b = mathieu_b(m,q)
    fc = specfun.fcoef(kd,m,q,b)
    return fc[:km]


def lpmn(m,n,z):
    """Associated Legendre functions of the first kind, Pmn(z) and its
    derivative, ``Pmn'(z)`` of order m and degree n.  Returns two
    arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and ``Pmn'(z)`` for
    all orders from ``0..m`` and degrees from ``0..n``.

    Parameters
    ----------
    m : int
       ``|m| <= n``; the order of the Legendre function.
    n : int
       where ``n >= 0``; the degree of the Legendre function.  Often
       called ``l`` (lower case L) in descriptions of the associated
       Legendre function
    z : float or complex
        Input value.

    Returns
    -------
    Pmn_z : (m+1, n+1) array
       Values for all orders 0..m and degrees 0..n
    Pmn_d_z : (m+1, n+1) array
       Derivatives for all orders 0..m and degrees 0..n
    """
    if not isscalar(m) or (abs(m)>n):
        raise ValueError("m must be <= n.")
    if not isscalar(n) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if not isscalar(z):
        raise ValueError("z must be scalar.")
    if (m < 0):
        mp = -m
        mf,nf = mgrid[0:mp+1,0:n+1]
        sv = errprint(0)
        fixarr = where(mf>nf,0.0,(-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
        sv = errprint(sv)
    else:
        mp = m
    if iscomplex(z):
        p,pd = specfun.clpmn(mp,n,real(z),imag(z))
    else:
        p,pd = specfun.lpmn(mp,n,z)
    if (m < 0):
        p = p * fixarr
        pd = pd * fixarr
    return p,pd


def lqmn(m,n,z):
    """Associated Legendre functions of the second kind, Qmn(z) and its
    derivative, ``Qmn'(z)`` of order m and degree n.  Returns two
    arrays of size ``(m+1, n+1)`` containing ``Qmn(z)`` and ``Qmn'(z)`` for
    all orders from ``0..m`` and degrees from ``0..n``.

    z can be complex.
    """
    if not isscalar(m) or (m<0):
        raise ValueError("m must be a non-negative integer.")
    if not isscalar(n) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if not isscalar(z):
        raise ValueError("z must be scalar.")
    m = int(m)
    n = int(n)

    # Ensure neither m nor n == 0
    mm = max(1,m)
    nn = max(1,n)

    if iscomplex(z):
        q,qd = specfun.clqmn(mm,nn,z)
    else:
        q,qd = specfun.lqmn(mm,nn,z)
    return q[:(m+1),:(n+1)],qd[:(m+1),:(n+1)]


def bernoulli(n):
    """Return an array of the Bernoulli numbers B0..Bn
    """
    if not isscalar(n) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    n = int(n)
    if (n < 2): n1 = 2
    else: n1 = n
    return specfun.bernob(int(n1))[:(n+1)]

def euler(n):
    """Return an array of the Euler numbers E0..En (inclusive)
    """
    if not isscalar(n) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    n = int(n)
    if (n < 2): n1 = 2
    else:  n1 = n
    return specfun.eulerb(n1)[:(n+1)]

def lpn(n,z):
    """Compute sequence of Legendre functions of the first kind (polynomials),
    Pn(z) and derivatives for all degrees from 0 to n (inclusive).

    See also special.legendre  for polynomial class.
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n < 1): n1 = 1
    else: n1 = n
    if iscomplex(z):
        pn,pd = specfun.clpn(n1,z)
    else:
        pn,pd = specfun.lpn(n1,z)
    return pn[:(n+1)],pd[:(n+1)]

## lpni

def lqn(n,z):
    """Compute sequence of Legendre functions of the second kind,
    Qn(z) and derivatives for all degrees from 0 to n (inclusive).
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (n!= floor(n)) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if (n < 1): n1 = 1
    else: n1 = n
    if iscomplex(z):
        qn,qd = specfun.clqn(n1,z)
    else:
        qn,qd = specfun.lqnb(n1,z)
    return qn[:(n+1)],qd[:(n+1)]

def ai_zeros(nt):
    """Compute the zeros of Airy Functions Ai(x) and Ai'(x), a and a'
    respectively, and the associated values of Ai(a') and Ai'(a).

    Returns
    -------
    a[l-1]   -- the lth zero of Ai(x)
    ap[l-1]  -- the lth zero of Ai'(x)
    ai[l-1]  -- Ai(ap[l-1])
    aip[l-1] -- Ai'(a[l-1])
    """
    kf = 1
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be a positive integer scalar.")
    return specfun.airyzo(nt,kf)

def bi_zeros(nt):
    """Compute the zeros of Airy Functions Bi(x) and Bi'(x), b and b'
    respectively, and the associated values of Ai(b') and Ai'(b).

    Returns
    -------
    b[l-1]   -- the lth zero of Bi(x)
    bp[l-1]  -- the lth zero of Bi'(x)
    bi[l-1]  -- Bi(bp[l-1])
    bip[l-1] -- Bi'(b[l-1])
    """
    kf = 2
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be a positive integer scalar.")
    return specfun.airyzo(nt,kf)

def lmbda(v,x):
    """Compute sequence of lambda functions with arbitrary order v
    and their derivatives.  Lv0(x)..Lv(x) are computed with v0=v-int(v).
    """
    if not (isscalar(v) and isscalar(x)):
        raise ValueError("arguments must be scalars.")
    if (v<0):
        raise ValueError("argument must be > 0.")
    n = int(v)
    v0 = v - n
    if (n < 1): n1 = 1
    else: n1 = n
    v1 = n1 + v0
    if (v!=floor(v)):
        vm, vl, dl = specfun.lamv(v1,x)
    else:
        vm, vl, dl = specfun.lamn(v1,x)
    return vl[:(n+1)], dl[:(n+1)]

def pbdv_seq(v,x):
    """Compute sequence of parabolic cylinder functions Dv(x) and
    their derivatives for Dv0(x)..Dv(x) with v0=v-int(v).
    """
    if not (isscalar(v) and isscalar(x)):
        raise ValueError("arguments must be scalars.")
    n = int(v)
    v0 = v-n
    if (n < 1): n1=1
    else: n1 = n
    v1 = n1 + v0
    dv,dp,pdf,pdd = specfun.pbdv(v1,x)
    return dv[:n1+1],dp[:n1+1]

def pbvv_seq(v,x):
    """Compute sequence of parabolic cylinder functions Dv(x) and
    their derivatives for Dv0(x)..Dv(x) with v0=v-int(v).
    """
    if not (isscalar(v) and isscalar(x)):
        raise ValueError("arguments must be scalars.")
    n = int(v)
    v0 = v-n
    if (n <= 1): n1=1
    else: n1 = n
    v1 = n1 + v0
    dv,dp,pdf,pdd = specfun.pbvv(v1,x)
    return dv[:n1+1],dp[:n1+1]

def pbdn_seq(n,z):
    """Compute sequence of parabolic cylinder functions Dn(z) and
    their derivatives for D0(z)..Dn(z).
    """
    if not (isscalar(n) and isscalar(z)):
        raise ValueError("arguments must be scalars.")
    if (floor(n)!=n):
        raise ValueError("n must be an integer.")
    if (abs(n) <= 1):
        n1 = 1
    else:
        n1 = n
    cpb,cpd = specfun.cpbdn(n1,z)
    return cpb[:n1+1],cpd[:n1+1]

def ber_zeros(nt):
    """Compute nt zeros of the kelvin function ber x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,1)

def bei_zeros(nt):
    """Compute nt zeros of the kelvin function bei x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,2)

def ker_zeros(nt):
    """Compute nt zeros of the kelvin function ker x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,3)

def kei_zeros(nt):
    """Compute nt zeros of the kelvin function kei x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,4)

def berp_zeros(nt):
    """Compute nt zeros of the kelvin function ber' x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,5)

def beip_zeros(nt):
    """Compute nt zeros of the kelvin function bei' x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,6)

def kerp_zeros(nt):
    """Compute nt zeros of the kelvin function ker' x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,7)

def keip_zeros(nt):
    """Compute nt zeros of the kelvin function kei' x
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,8)

def kelvin_zeros(nt):
    """Compute nt zeros of all the kelvin functions returned in a
    length 8 tuple of arrays of length nt.
    The tuple containse the arrays of zeros of
    (ber, bei, ker, kei, ber', bei', ker', kei')
    """
    if not isscalar(nt) or (floor(nt)!=nt) or (nt<=0):
        raise ValueError("nt must be positive integer scalar.")
    return specfun.klvnzo(nt,1), \
           specfun.klvnzo(nt,2), \
           specfun.klvnzo(nt,3), \
           specfun.klvnzo(nt,4), \
           specfun.klvnzo(nt,5), \
           specfun.klvnzo(nt,6), \
           specfun.klvnzo(nt,7), \
           specfun.klvnzo(nt,8)

def pro_cv_seq(m,n,c):
    """Compute a sequence of characteristic values for the prolate
    spheroidal wave functions for mode m and n'=m..n and spheroidal
    parameter c.
    """
    if not (isscalar(m) and isscalar(n) and isscalar(c)):
        raise ValueError("Arguments must be scalars.")
    if (n!=floor(n)) or (m!=floor(m)):
        raise ValueError("Modes must be integers.")
    if (n-m > 199):
        raise ValueError("Difference between n and m is too large.")
    maxL = n-m+1
    return specfun.segv(m,n,c,1)[1][:maxL]

def obl_cv_seq(m,n,c):
    """Compute a sequence of characteristic values for the oblate
    spheroidal wave functions for mode m and n'=m..n and spheroidal
    parameter c.
    """
    if not (isscalar(m) and isscalar(n) and isscalar(c)):
        raise ValueError("Arguments must be scalars.")
    if (n!=floor(n)) or (m!=floor(m)):
        raise ValueError("Modes must be integers.")
    if (n-m > 199):
        raise ValueError("Difference between n and m is too large.")
    maxL = n-m+1
    return specfun.segv(m,n,c,-1)[1][:maxL]

def ellipk(m):
    """y=ellipk(m) returns the complete integral of the first kind:
    integral(1/sqrt(1-m*sin(t)**2),t=0..pi/2)

    This function is rather imprecise around m==1. For more precision
    around this point, use ellipkm1."""
    return ellipkm1(1 - asarray(m))

def agm(a,b):
    """Arithmetic, Geometric Mean

    Start with a_0=a and b_0=b and iteratively compute

    a_{n+1} = (a_n+b_n)/2
    b_{n+1} = sqrt(a_n*b_n)

    until a_n=b_n.   The result is agm(a,b)

    agm(a,b)=agm(b,a)
    agm(a,a) = a
    min(a,b) < agm(a,b) < max(a,b)
    """
    s = a + b + 0.0
    return (pi / 4) * s / ellipkm1(4 * a * b / s ** 2)