QM: Hydrogen radial wavefunctions added

Tests and doctests added. The wavefunctions were carefully checked one by one
for all (n, l) with n <= 4 and all those are tested. We also test the
normalization for all n <= 2.

sympy/physics/hydrogen.py

@@ -0,0 +1,52 @@
+from sympy import factorial, sqrt, exp, S, laguerre_l
+
+def R_nl(n, l, a, r):
+    """
+    Returns the Hydrogen radial wavefunction R_{nl}.
+
+    n, l .... quantum numbers 'n' and 'l'
+    a    .... Bohr radius \hbar^2 / (Z*M)
+    r    .... radial coordinate
+
+    Examples::
+
+    >>> from sympy.physics.hydrogen import R_nl
+    >>> from sympy import var
+    >>> var("r a")
+    (r, a)
+    >>> R_nl(1, 0, a, r)
+    2*(a**(-3))**(1/2)*exp(-r/a)
+    >>> R_nl(2, 0, a, r)
+    2**(1/2)*(a**(-3))**(1/2)*(2 - r/a)*exp(-r/(2*a))/4
+    >>> R_nl(2, 1, a, r)
+    r*6**(1/2)*(a**(-3))**(1/2)*exp(-r/(2*a))/(12*a)
+
+    In most cases you probably want to set a=1::
+
+    >>> R_nl(1, 0, 1, r)
+    2*exp(-r)
+    >>> R_nl(2, 0, 1, r)
+    2**(1/2)*(2 - r)*exp(-r/2)/4
+    >>> R_nl(3, 0, 1, r)
+    2*3**(1/2)*(3 - 2*r + 2*r**2/9)*exp(-r/3)/27
+
+    The normalization of the radial wavefunction is::
+
+    >>> from sympy import integrate, oo
+    >>> integrate(R_nl(1, 0, 1, r)**2 * r**2, (r, 0, oo))
+    1
+    >>> integrate(R_nl(2, 0, 1, r)**2 * r**2, (r, 0, oo))
+    1
+    >>> integrate(R_nl(2, 1, 1, r)**2 * r**2, (r, 0, oo))
+    1
+
+    """
+    # sympify arguments
+    n, l, a, r = S(n), S(l), S(a), S(r)
+    # radial quantum number
+    n_r = n - l - 1
+    # rescaled "r"
+    r0 = 2 * r / (n * a)
+    # normalization coefficient
+    C =  sqrt((S(2)/(n*a))**3 * factorial(n_r) / (2*n*factorial(n+l)))
+    return  C * r0**l * laguerre_l(n_r, 2*l+1, r0) * exp(-r0/2)

sympy/physics/tests/test_hydrogen.py

@@ -0,0 +1,34 @@
+from sympy import var, sqrt, exp, simplify, S, integrate, oo
+from sympy.physics.hydrogen import R_nl
+
+var("a r")
+
+def test_wavefunction():
+    R = {
+            (1, 0): 2*sqrt(1/a**3) * exp(-r/a),
+            (2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1-r/(2*a)),
+            (2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
+            (3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) * \
+                    (1-2*r/(3*a) + S(2)/27 * (r/a)**2),
+            (3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) * \
+                    (1-r/(6*a)) * r/a,
+            (3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
+            (4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) * \
+                    (1-3*r/(4*a)+S(1)/8 * (r/a)**2-S(1)/192 * (r/a)**3),
+            (4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) * \
+                    (1-r/(4*a)+S(1)/80 * (r/a)**2) * (r/a),
+            (4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) * \
+                    (1-r/(12*a)) * (r/a)**2,
+            (4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
+            }
+    for n, l in R:
+        assert simplify(R_nl(n, l, a, r) - R[(n, l)]) == 0
+
+def test_norm():
+    # Maximum "n" which is tested:
+    n_max = 2
+    # you can test any n and it works, but it's slow, so it's commented out:
+    #n_max = 4
+    for n in range(n_max+1):
+        for l in range(n):
+            assert integrate(R_nl(n, l, 1, r)**2 * r**2, (r, 0, oo)) == 1