"Spontaneous emission coefficients for H-alpha"

psi(n,l,m) = R(n,l) Y(l,m)

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))

-- associated Laguerre polynomial (k is a local var)

L(x,n,m,k) = (n + m)! sum(k, 0, n, (-x)^k / ((n - k)! (m + k)! k!))

-- spherical harmonic

Y(l,m) = (-1)^m sqrt((2 l + 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),
           (sin(theta)/2)^m sum(k, 0, l - m,
           (-1)^k (l + m + k)! / (l - m - k)! / (m + k)! / k! *
           ((1 - cos(theta)) / 2)^k))

-- integrate f

I(f) = do(
  f = f r^2 sin(theta), -- multiply by volume element
  f = expform(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))

-- average over 3 initial states per final state

X3p2s = 1/3 X3p2s

-- average over 5 initial states per final state

X3d2p = 1/5 X3d2p

E(n) = -alpha hbar c / (2 n^2 a0)
omega32 = (E(3) - E(2)) / hbar

A3s2p = 4 alpha omega32^3 X3s2p / (3 c^2)
A3p2s = 4 alpha omega32^3 X3p2s / (3 c^2)
A3d2p = 4 alpha omega32^3 X3d2p / (3 c^2)

-- CODATA Internationally recommended 2022 values
-- https://physics.nist.gov/cuu/Constants/
-- c, e, h, and k are exact values

a0 = 5.29177210544 10^(-11) meter
alpha = 7.2973525643 10^(-3)
c = 299792458.0 meter / second
e = 1.602176634 10^(-19) coulomb
epsilon0 = 8.8541878188 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.1093837139 10^(-31) kilogram
mp = 1.67262192595 10^(-27) kilogram
mu0 = 1.25663706127 10^(-6) newton / ampere^2

coulomb = ampere second
farad = coulomb / volt
joule = kilogram meter^2 / second^2
newton = kilogram meter / second^2
tesla = kilogram / second^2 / ampere
volt = joule / coulomb

ampere = "ampere"
kelvin = "kelvin"
kilogram = "kilogram"
meter = "meter"
second = "second"

pi = float(pi) -- use numerical value of pi
mu = me mp / (me + mp)
a0 = a0 me / mu -- correction for reduced electron mass

A3s2p
A3p2s
A3d2p

"Verify coefficients"

err(a,b) = 2 abs((a - b) / (a + b)) -- relative error

check(err(A3s2p, 6.313 10^6 / second) < 0.0001)
check(err(A3p2s, 2.245 10^7 / second) < 0.0001)
check(err(A3d2p, 6.465 10^7 / second) < 0.0001)

"ok"