sho.py 2.4 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495
  1. from sympy.core import S, pi, Rational
  2. from sympy.functions import assoc_laguerre, sqrt, exp, factorial, factorial2
  3. def R_nl(n, l, nu, r):
  4. """
  5. Returns the radial wavefunction R_{nl} for a 3d isotropic harmonic
  6. oscillator.
  7. Parameters
  8. ==========
  9. ``n`` :
  10. The "nodal" quantum number. Corresponds to the number of nodes in
  11. the wavefunction. ``n >= 0``
  12. ``l`` :
  13. The quantum number for orbital angular momentum.
  14. ``nu`` :
  15. mass-scaled frequency: nu = m*omega/(2*hbar) where `m` is the mass
  16. and `omega` the frequency of the oscillator.
  17. (in atomic units ``nu == omega/2``)
  18. ``r`` :
  19. Radial coordinate.
  20. Examples
  21. ========
  22. >>> from sympy.physics.sho import R_nl
  23. >>> from sympy.abc import r, nu, l
  24. >>> R_nl(0, 0, 1, r)
  25. 2*2**(3/4)*exp(-r**2)/pi**(1/4)
  26. >>> R_nl(1, 0, 1, r)
  27. 4*2**(1/4)*sqrt(3)*(3/2 - 2*r**2)*exp(-r**2)/(3*pi**(1/4))
  28. l, nu and r may be symbolic:
  29. >>> R_nl(0, 0, nu, r)
  30. 2*2**(3/4)*sqrt(nu**(3/2))*exp(-nu*r**2)/pi**(1/4)
  31. >>> R_nl(0, l, 1, r)
  32. r**l*sqrt(2**(l + 3/2)*2**(l + 2)/factorial2(2*l + 1))*exp(-r**2)/pi**(1/4)
  33. The normalization of the radial wavefunction is:
  34. >>> from sympy import Integral, oo
  35. >>> Integral(R_nl(0, 0, 1, r)**2*r**2, (r, 0, oo)).n()
  36. 1.00000000000000
  37. >>> Integral(R_nl(1, 0, 1, r)**2*r**2, (r, 0, oo)).n()
  38. 1.00000000000000
  39. >>> Integral(R_nl(1, 1, 1, r)**2*r**2, (r, 0, oo)).n()
  40. 1.00000000000000
  41. """
  42. n, l, nu, r = map(S, [n, l, nu, r])
  43. # formula uses n >= 1 (instead of nodal n >= 0)
  44. n = n + 1
  45. C = sqrt(
  46. ((2*nu)**(l + Rational(3, 2))*2**(n + l + 1)*factorial(n - 1))/
  47. (sqrt(pi)*(factorial2(2*n + 2*l - 1)))
  48. )
  49. return C*r**(l)*exp(-nu*r**2)*assoc_laguerre(n - 1, l + S.Half, 2*nu*r**2)
  50. def E_nl(n, l, hw):
  51. """
  52. Returns the Energy of an isotropic harmonic oscillator.
  53. Parameters
  54. ==========
  55. ``n`` :
  56. The "nodal" quantum number.
  57. ``l`` :
  58. The orbital angular momentum.
  59. ``hw`` :
  60. The harmonic oscillator parameter.
  61. Notes
  62. =====
  63. The unit of the returned value matches the unit of hw, since the energy is
  64. calculated as:
  65. E_nl = (2*n + l + 3/2)*hw
  66. Examples
  67. ========
  68. >>> from sympy.physics.sho import E_nl
  69. >>> from sympy import symbols
  70. >>> x, y, z = symbols('x, y, z')
  71. >>> E_nl(x, y, z)
  72. z*(2*x + y + 3/2)
  73. """
  74. return (2*n + l + Rational(3, 2))*hw