## Automatically adapted for scipy Oct 05, 2005 by convertcode.py

#
# Author:  Travis Oliphant, 2002
#

from numpy import *
from _cephes import *
import types
import specfun

def sinc(x):
    """Returns sin(pi*x)/(pi*x) at all points of array x.
    """
    w = asarray(asarray(x)*pi)
    return where(x==0, 1.0, sin(w)/w)

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.

    Outputs (all are arrays of length nt):

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

    See 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[:nt],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 >= 0
        The degree of the harmonic.
      - `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`.
    """
    x = cos(phi)
    m,n = int(m), int(n)
    Pmn,Pmnd = lpmn(m,n,x)
    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):
    gam = gamma
    fac = gam(k+1+n)/gam(k+1)/gam(n+1)
    return fac*hyp1f1(-n,k+1,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.

    z can be complex.
    """
    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).

    Outputs:

      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).

    Outputs:

      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 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)
    """
    res1 = a+b+0.0
    res2 = a-b
    k = res2 / res1
    return res1*pi/4/ellipk(k**2)