The monopoles use an expression like hn(k*rs) * bn(k*r), which exist as a limit k->0 but is undefined for k==0. Consequently spherical_ps and circular_ls yield NaN for k==0.
The DC case must be handled separately.
Hopefully, a closed-form expression for the limit can be found.