assume(a > 0, e > 0, ε_0 > 0, n > 0, l > 0);

assume(0 < φ, φ < 2*π, 0 < θ, θ < π, r > 0);

fpprintprec : 4;

simp(e) := ratsubst(abs(sin(θ)), sqrt(1-cos(θ))*sqrt(cos(θ)+1), trigsimp(radcan(ratsimp(e))));

ndiff(expr, n) := (
    (α : expr), thru n do α: simp(diff(α, x)), α);

Pl(l) := simp(1/(2^l * l!) * ndiff( (x^2 - 1)^l, l));

Pml(m, l) := simp((-1)^m * (1 - x^2)^(abs(m) / 2) * ndiff(Pl(l), abs(m)));

Yml(m, l) := simp((-1)^m * sqrt( (2 * l + 1) / (4*π) * (l - abs(m))! / (l + abs(m))! ) * %e^(%i *m * φ) * subst(cos(θ), x, Pml(m, l)));

Lq(q) := simp(%e^x * ndiff(%e^-x * x^q, q));

Lpq(p, q) := simp((-1)^p * ndiff(Lq(q + p), p));

norm_factor(n, l, m) :=  simp(sqrt( (2 / (n*a))^3 * (n - l - 1)! / ( 2* n * ((n+l)!)^3 ) ) * %e^(-r / (n*a)) * (2*r / (n*a))^l);

ψ(n, l, m) := simp(norm_factor(n, l, m) * subst(2*r/(n*a), x, Lpq(2*l + 1, n - l - 1)) * Yml(m, l));

rv : [ r * sin(θ) * cos(φ), r * sin(θ) * sin(φ), r*cos(θ)];

spherical_int(expr) := integrate(integrate(integrate(expr * r^2*sin(θ), r, 0, inf), φ, 0, 2*π), θ, 0, π);

ρ(n1, m1, l1, n2, m2, l2) := e * spherical_int(conjugate(ψ(n2, m2, l2)) * rv * ψ(n1, m1, l1));

norm(x) := sqrt(x . conjugate(x));

En(n) := -(m / (2*ħ^2)) * (e^2 / (4*π*ε_0))^2 / n^2;

ω_0(n1, n2) := simp((En(n1) - En(n2))  / ħ);

A(n1, l1, m1, n2, m2, l2) := simp((ω_0(n1, n2))^3 * norm(ρ(n1, l1, m1, n2, m2, l2))^2 /  (3*π * ε_0 * ħ * c^3));

Af(n1, l1, m1, n2, m2, l2) := float(sublis(bfloat([c=%c, ħ=%ħ, ε_0=%ε_0, a=%a, e=%el, m=%m_e]), A(n1, l1, m1, n2, m2, l2)));

l2_values(n1, l1, n2) := sublist([l1 - 1, l1 + 1], lambda([v], is(v >= 0 and v < n2)));

m2_values(n1, m1, n2, l2) := sublist([m1 - 1, m1, m1 + 1], lambda([v], is(v >= -l2 and v <= l2)));

m2_values(4, 2, 3, 2);

decay_modes_with_trans(n1, l1, m1) :=
    (dm: [],
    for n2 from 1 thru n1-1 do
    for l2 in l2_values(n1, l1, n2) do
    for m2 in m2_values(n1, m1, n2, l2) do
    push( [n2, l2, m2, Af(n1, l1, m1, n2, l2, m2)], dm),
    dm)
        ;

lifetime_from_modes(modes) := (s: lsum(m[4], m, modes), if is(bfloat(s) = 0b0) then inf else 1/s);

texfloat(f) := concat("\\(", tex1(f), "\\)");

assume(sin(θ) > 0);

for n1 from 1 thru 6 do
for l1 from 0 thru n1-1 do
for m1 from -l1 thru l1 do
(print("<tr><td>"),
print(concat(n1, ", ", l1, ", ", m1), "</td><td><div>"),    
print("\\(", tex1(ψ(n1, l1, m1)), "\\)"),
print("</div></td><td>"),
dm : decay_modes_with_trans(n1, l1, m1),
for mode in dm do
print(concat(mode[1], ",", mode[2], ",", mode[3]), "</br>"),
print("</td><td>"),
for mode in dm do
print(texfloat(mode[4]), "<br/>"),
print("</td><td>"),
print(texfloat(lifetime_from_modes(dm))),
print("</td></tr>"));