def wigner_d_naive_v3(l, m, n, approx_lim=1000000):
    """
    Wigner "small d" matrix. (Euler z-y-z convention)
    example:
        l = 2
        m = 1
        n = 0
        beta = linspace(0,pi,100)
        wd210 = wignerd(l,m,n)(beta)

    some conditions have to be met:
         l >= 0
        -l <= m <= l
        -l <= n <= l

    The approx_lim determines at what point
    bessel functions are used. Default is when:
        l > m+10
          and
        l > n+10

    for integer l and n=0, we can use the spherical harmonics. If in
    addition m=0, we can use the ordinary legendre polynomials.
    """
    from scipy.special import jv, legendre, sph_harm, jacobi
    try:
        from scipy.misc import factorial, comb
    except:
        from scipy.special import factorial, comb
    from numpy import floor, sqrt, sin, cos, exp, power
    from math import pi
    from scipy.special import jacobi
    if l < 0 or abs(m) > l or abs(n) > l:
        raise ValueError('wignerd(l = {0}, m = {1}, n = {2}) value error.'.
            format(l, m, n) + ' Valid range for parameters: l>=0, -l<=m,n<=l.')
    if l > m + approx_lim and l > n + approx_lim:
        return lambda beta: jv(m - n, l * beta)
    if floor(l) == l and n == 0:
        if m == 0:
            return lambda beta: legendre(l)(cos(beta))
        elif False:
            a = sqrt(4.0 * pi / (2.0 * l + 1.0))
            return lambda beta: a * sph_harm(m, l, beta, 0.0).conj()
    jmn_terms = {(l + n): (m - n, m - n), (l - n): (n - m, 0.0), (l + m): (
        n - m, 0.0), (l - m): (m - n, m - n)}
    k = min(jmn_terms)
    a, lmb = jmn_terms[k]
    b = 2.0 * l - 2.0 * k - a
    if a < 0 or b < 0:
        raise ValueError('wignerd(l = {0}, m = {1}, n = {2}) value error.'.
            format(l, m, n) +
            ' Encountered negative values in (a,b) = ({0},{1})'.format(a, b))
    coeff = power(-1.0, lmb) * sqrt(comb(2.0 * l - k, k + a)) * (1.0 / sqrt
        (comb(k + b, b)))
    return lambda beta: coeff * power(sin(0.5 * beta), a) * power(cos(0.5 *
        beta), b) * jacobi(k, a, b)(cos(beta))
