[Bug] There is a small bug about function d_me

[floatingCatty]

The d_me fuction:

def d_me(N, l, m1, m2): """Computes rotational matrix elements according to A.V. Podolskiy and P. Vogl, Phys. Rev. B. 69, 233101 (2004) Parameters

N :
    directional cosine relative to z-axis
l :
    orbital quantum number
m1 :
    magnetic quantum number
m2 :
    magnetic quantum number
Returns
-------
type
    rotational matrix element
"""

if N == -1.0 and m1 == m2:
    prefactor = math.sqrt(math.factorial(l + m2) * math.factorial(l - m2) *
                          math.factorial(l + m1) * math.factorial(l - m1))
else:
    prefactor = ((0.5 * (1 + N)) ** l) * (((1 - N) / (1 + N)) ** (m1 * 0.5 - m2 * 0.5)) * \
        math.sqrt(math.factorial(l + m2) * math.factorial(l - m2) *
                  math.factorial(l + m1) * math.factorial(l - m1))

ans = 0
for t in range(2 * l + 2):
    if l + m2 - t >= 0 and l - m1 - t >= 0 and t + m1 - m2 >= 0:
        if N == -1.0 and t == 0:
            ans += ((-1) ** t) / \
                   (math.factorial(l + m2 - t) * math.factorial(l - m1 - t) *
                    math.factorial(t) * math.factorial(t + m1 - m2))
        else:
            ans += ((-1) ** t) * (((1 - N) / (1 + N)) ** t) / \
                   (math.factorial(l + m2 - t) * math.factorial(l - m1 - t) *
                    math.factorial(t) * math.factorial(t + m1 - m2))

return np.nan_to_num(ans * prefactor)

if N, l, m1, m2 equals -1, 1, 0, 0 respectively, then (1 - N) / (1 + N) in the "and +=" line will raise an error of dividing by zero. should there be a condition that excludes the N = -1?

[freude]

Thank you very much for pointing this out. I have fixed this by adding an infinitesimal constant to N when N=-1. This fix has passed all the tests.