-- psi returns a hydrogen atom eigenfunction psi(n,l,m) = R(n,l) Y(l,m) -- R returns a radial eigenfunction R(n,l) = 2 / n^2 * a0^(-3/2) * sqrt((n - l - 1)! / (n + l)!) * (2 r / (n a0))^l * L(2 r / (n a0),n - l - 1,2 l + 1) * exp(-r / (n a0)) L(x,n,m) = (n + m)! sum(k,0,n,(-x)^k / ((n - k)! (m + k)! k!)) -- Y returns a spherical harmonic eigenfunction Y(l,m) = (-1)^m sqrt((2l + 1) / (4 pi) (l - m)! / (l + m)!) * P(cos(theta),l,m) exp(i m phi) P(f,n,m) = eval(1 / (2^n n!) (1 - x^2)^(m/2) * d((x^2 - 1)^n,x,n + m),x,f) -- E(n) returns the nth energy eigenvalue E(n) = -e^2 / (8 pi epsilon0 a0 n^2) -- integrate f I(f) = do( f = eval(f,sqrt(1 - cos(theta)^2),sin(theta)), -- simplify f = eval(f,1/sqrt(1 - cos(theta)^2),1/sin(theta)), -- simplify f = f r^2 sin(theta), -- multiply by volume element f = defint(f,theta,0,pi,phi,0,2pi), f = integral(f,r), 0 - eval(f,r,0) ) xbar1 = I(conj(psi(1,0,0)) r sin(theta) cos(phi) psi(2,1,1)) xbar2 = I(conj(psi(1,0,0)) r sin(theta) cos(phi) psi(2,1,-1)) xbar3 = I(conj(psi(1,0,0)) r sin(theta) cos(phi) psi(2,1,0)) xbar4 = I(conj(psi(1,0,0)) r sin(theta) cos(phi) psi(2,0,0)) ybar1 = I(conj(psi(1,0,0)) r sin(theta) sin(phi) psi(2,1,1)) ybar2 = I(conj(psi(1,0,0)) r sin(theta) sin(phi) psi(2,1,-1)) ybar3 = I(conj(psi(1,0,0)) r sin(theta) sin(phi) psi(2,1,0)) ybar4 = I(conj(psi(1,0,0)) r sin(theta) sin(phi) psi(2,0,0)) zbar1 = I(conj(psi(1,0,0)) r cos(theta) psi(2,1,1)) zbar2 = I(conj(psi(1,0,0)) r cos(theta) psi(2,1,-1)) zbar3 = I(conj(psi(1,0,0)) r cos(theta) psi(2,1,0)) zbar4 = I(conj(psi(1,0,0)) r cos(theta) psi(2,0,0)) rbar1 = conj(xbar1) xbar1 + conj(ybar1) ybar1 + conj(zbar1) zbar1 rbar2 = conj(xbar2) xbar2 + conj(ybar2) ybar2 + conj(zbar2) zbar2 rbar3 = conj(xbar3) xbar3 + conj(ybar3) ybar3 + conj(zbar3) zbar3 rbar4 = conj(xbar4) xbar4 + conj(ybar4) ybar4 + conj(zbar4) zbar4 -- print psi = quote(psi) (("",psi(2,1,1),psi(2,1,-1),psi(2,1,0),psi(2,0,0)), (x,xbar1,xbar2,xbar3,xbar4), (y,ybar1,ybar2,ybar3,ybar4), (z,zbar1,zbar2,zbar3,zbar4), (r^2,rbar1,rbar2,rbar3,rbar4)) -- physical constants (h and k are exact values) c = 299792458.0 meter / second e = 1.602176634 10^(-19) coulomb epsilon0 = 8.8541878128 10^(-12) farad / meter h = 6.62607015 10^(-34) joule second hbar = h / float(2 pi) k = 1.380649 10^(-23) joule / kelvin me = 9.1093837015 10^(-31) kilogram mp = 1.67262192369 10^(-27) kilogram mu = me mp / (me + mp) -- derived units coulomb = ampere second farad = coulomb / volt joule = kilogram meter^2 / second^2 volt = joule / coulomb -- base units (for printing) ampere = "ampere" kelvin = "kelvin" kilogram = "kilogram" meter = "meter" second = "second" pi = float(pi) a0 = 4 pi epsilon0 hbar^2 / (e^2 me) a0 omega21 = 1/hbar (E(2) - E(1)) omega21 A21 = e^10 me / (26244 pi^5 epsilon0^5 hbar^6 c^3) A21 B21 = c^2 / (2 h nu^3) A21 B21 g1 = 2 g2 = 6 B12 = g2 / g1 B21 B12