-- Compute matrices for 2d hydrogen E(n) = -k^2 mu / (2 hbar^2 (n + 1/2)^2) e = quote(e) -- clear e k = e^2 / (4 pi epsilon0) H = ((E(0),0,0,0), (0,E(1),0,0), (0,0,E(2),0), (0,0,0,E(3))) psi(n,m) = 1 / a0 / sqrt(pi) * (n + 1/2)^(-3/2) * sqrt((n - abs(m))! / (n + abs(m))!) * (2 r / a0 / (n + 1/2))^abs(m) * Laguerre(2 r / a0 / (n + 1/2), n - abs(m), 2 abs(m)) * exp(-r / a0 / (n + 1/2)) * exp(i m phi) a0 = hbar^2 / (k mu) -- Laguerre polynomial (y is a local variable) Laguerre(x,n,m,y) = eval(y^(-m) exp(y) / n! d(exp(-y) y^(n + m), y, n), y, x) psi0 = psi(0,0) psi1 = psi(1,0) psi2 = psi(2,0) psi3 = psi(3,0) I(f) = do( f = defint(f r, phi, 0, 2 pi), f = integral(f,r), 0 - eval(f,r,0) ) phat(psi) = -hbar^2 (d(psi,r,2) + d(psi,r) / r + d(psi,phi,2) / r^2) P2 = ((I(conj(psi0) phat(psi0)), I(conj(psi0) phat(psi1)), I(conj(psi0) phat(psi2)), I(conj(psi0) phat(psi3))), (I(conj(psi1) phat(psi0)), I(conj(psi1) phat(psi1)), I(conj(psi1) phat(psi2)), I(conj(psi1) phat(psi3))), (I(conj(psi2) phat(psi0)), I(conj(psi2) phat(psi1)), I(conj(psi2) phat(psi2)), I(conj(psi2) phat(psi3))), (I(conj(psi3) phat(psi0)), I(conj(psi3) phat(psi1)), I(conj(psi3) phat(psi2)), I(conj(psi3) phat(psi3)))) V = ((I(conj(psi0) 1/r psi0), I(conj(psi0) 1/r psi1), I(conj(psi0) 1/r psi2), I(conj(psi0) 1/r psi3)), (I(conj(psi1) 1/r psi0), I(conj(psi1) 1/r psi1), I(conj(psi1) 1/r psi2), I(conj(psi1) 1/r psi3)), (I(conj(psi2) 1/r psi0), I(conj(psi2) 1/r psi1), I(conj(psi2) 1/r psi2), I(conj(psi2) 1/r psi3)), (I(conj(psi3) 1/r psi0), I(conj(psi3) 1/r psi1), I(conj(psi3) 1/r psi2), I(conj(psi3) 1/r psi3))) P2 = simplify(P2) V = simplify(V) H P2 V "Verify symbolic matrices (1=ok)" H == P2 / (2 mu) - k V "Compute numerical values..." -- physical constants (h and kB 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) kB = 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) -- use numerical value of pi H P2 V