Repository: sympy/sympy
File: sympy/physics/hydrogen.py (Python)

Code before:

    >>> import sympy
    >>> from sympy.physics.hydrogen import Y_lm
    >>> phi, theta = sympy.symbols('phi theta', real=True)
    >>> Y_lm(1, 0, phi, theta)
    sqrt(3)*cos(theta)/(2*sqrt(pi))
    >>> Y_lm(2, 0, phi, theta).simplify()
    sqrt(5)*(3*cos(theta)**2 - 1)/(4*sqrt(pi))
    >>> Y_lm(1, -1, phi, theta, real=True).simplify()
    -sqrt(3)*sin(phi)*sin(theta)/(2*sqrt(pi))
    >>> Y_lm(1, -1, phi, theta, real=False).simplify()
    sqrt(6)*exp(-I*phi)*sin(theta)/(4*sqrt(pi))
    >>> Y_lm(1, 1, 1, 1).evalf()
    -0.157078470548806 - 0.244635223409689*I

    Find angular nodes in dz2 atomic orbital by solving for zero probability.
    This returns the node theta angles in radians.

    >>> dz2 = Y_lm(2, 0, phi, theta)**2
    >>> [sol.evalf() for sol in sympy.solve(dz2)]
    [4.0969092717143, 5.32786868905508, 2.18627603546528, 0.955316618124509]

    """
    l, m, theta, phi = map(S, [l, m, theta, phi])
    if m.is_integer and l.is_integer and abs(m) > l:
        raise ValueError("|'m'| must be less or equal 'l'")
    if real:
        return Znm(l, m, theta, phi).expand(func=True)
    else:
        return Ynm(l, m, theta, phi).expand(func=True)


def Psi_nlm(n, l, m, r, phi, theta, Z=1):
    """
    Returns the Hydrogen wave function psi_{nlm}. It's the product of
    the radial wavefunction R_{nl} and the spherical harmonic Y_{l}^{m}.

    Parameters
    ==========

    n : integer
        Principal Quantum Number which is
        an integer with possible values as 1, 2, 3, 4,...
    l : integer
        ``l`` is the Angular Momentum Quantum Number with
        values ranging from 0 to ``n-1``.
    m : integer
        ``m`` is the Magnetic Quantum Number with values
        ranging from ``-l`` to ``l``.
    r :
        radial coordinate

Review comment (bug):
This can still fail for symbols:
```python
In [7]: x = symbols('x', integer=True)

In [8]: x.is_integer
Out[8]: True

In [9]: if x > 1:
   ...:     print('greater')
   ...: 
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
...
TypeError: cannot determine truth value of Relational: x > 1
```
The usual way to write this would be:
```
if (abs(m) > l) is S.true:
```
That one can still fail for complex `l` but I don't know if that is relevant here.

Code after:
    sqrt(3)*cos(theta)/(2*sqrt(pi))
    >>> Y_lm(2, -1, phi, theta).simplify()
    sqrt(30)*exp(-I*phi)*sin(2*theta)/(8*sqrt(pi))
    >>> Y_lm(1, -1, phi, theta).simplify()
    sqrt(6)*exp(-I*phi)*sin(theta)/(4*sqrt(pi))
    >>> Y_lm(1, 1, 0, 1).evalf()
    -0.290723302201011

    Find angular nodes in dz2 atomic orbital by solving for zero probability.
    This returns the node theta angles in radians.

    >>> dz2 = Y_lm(2, 0, phi, theta)**2
    >>> [sol.evalf() for sol in sympy.solve(dz2)]
    [4.0969092717143, 5.32786868905508, 2.18627603546528, 0.955316618124509]

    See Also
    ========

    Z_lm, sympy.functions.special.spherical_harmonics.Ynm,
    sympy.functions.special.spherical_harmonics.Ynm_c,
    sympy.functions.special.spherical_harmonics.Znm

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Spherical_harmonics
    .. [2] https://en.wikipedia.org/wiki/Table_of_spherical_harmonics
    .. [3] https://mathworld.wolfram.com/SphericalHarmonic.html
    .. [4] https://functions.wolfram.com/Polynomials/SphericalHarmonicY/
    .. [5] https://dlmf.nist.gov/14.30

    """
    return Ynm(l, m, theta, phi).expand(func=True)


def Z_lm(l, m, phi, theta):
    """
    Returns the real spherical harmonic or angular function Y_{l}^{m}.

    Parameters
    ==========

    l : integer
        ``l`` is the Angular Momentum Quantum Number with
        values ranging from 0 to ``n-1``.
    m : iteger
        ``m`` is the Magnetic Quantum Number with values
        ranging from ``-l`` to ``l``.
    phi :
        azimuthal angle
    theta :