-- Compute spontaneous emission coefficients for H-alpha -- 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,k) = (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(l,m) exp(i m phi) -- associated Legendre of cos theta (arxiv.org/abs/1805.12125) P(l,m,k) = test(m < 0, (-1)^m (l + m)! / (l - m)! P(l,-m), 2^(-m) sin(theta)^m sum(k,0,l - m, (-1)^k (l + m + k)! / (l - m - k)! / (m + k)! / k! * ((1 - cos(theta)) / 2)^k)) -- E(n) returns the nth energy eigenvalue E(n) = -e^2 / (8 pi epsilon0 a0 n^2) -- integrate f I(f) = do( f = f r^2 sin(theta), -- multiply by volume element f = circexp(f), -- convert to exponential form f = defint(f,theta,0,pi,phi,0,2pi), f = integral(f,r), 0 - eval(f,r,0) -- return value ) X(fk,fi) = do( xki = I(conj(fk) r sin(theta) cos(phi) fi), yki = I(conj(fk) r sin(theta) sin(phi) fi), zki = I(conj(fk) r cos(theta) fi), conj(xki) xki + conj(yki) yki + conj(zki) zki -- return value ) X3s2p = X(psi(2,1,1),psi(3,0,0)) + X(psi(2,1,0),psi(3,0,0)) + X(psi(2,1,-1),psi(3,0,0)) X3p2s = X(psi(2,0,0),psi(3,1,1)) + X(psi(2,0,0),psi(3,1,0)) + X(psi(2,0,0),psi(3,1,-1)) X3d2p = X(psi(2,1,1),psi(3,2,2)) + X(psi(2,1,1),psi(3,2,1)) + X(psi(2,1,1),psi(3,2,0)) + X(psi(2,1,1),psi(3,2,-1)) + X(psi(2,1,1),psi(3,2,-1)) + X(psi(2,1,0),psi(3,2,2)) + X(psi(2,1,0),psi(3,2,1)) + X(psi(2,1,0),psi(3,2,0)) + X(psi(2,1,0),psi(3,2,-1)) + X(psi(2,1,0),psi(3,2,-1)) + X(psi(2,1,-1),psi(3,2,2)) + X(psi(2,1,-1),psi(3,2,1)) + X(psi(2,1,-1),psi(3,2,0)) + X(psi(2,1,-1),psi(3,2,-1)) + X(psi(2,1,-1),psi(3,2,-1)) -- physical constants (c, e, h, and k are exact values) c = 299792458.0 meter / second e = 1.602176634 10^(-19) coulomb h = 6.62607015 10^(-34) joule second hbar = h / float(2 pi) k = 1.380649 10^(-23) joule / kelvin epsilon0 = 8.8541878128 10^(-12) farad / meter 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 a0 = 4 pi epsilon0 hbar^2 / (e^2 mu) omega32 = (E(3) - E(2)) / hbar A3s2p = e^2 / (3 pi epsilon0 hbar c^3) omega32^3 X3s2p A3s2p A3p2s = e^2 / (3 pi epsilon0 hbar c^3) omega32^3 X3p2s / 3 A3p2s A3d2p = e^2 / (3 pi epsilon0 hbar c^3) omega32^3 X3d2p / 5 A3d2p Run