pub fn rho(s: Float, m1: Float, m2: Float) -> Complex<Float>
Computes the phase-space factor $\rho(s, m_1, m_2) = \sqrt(\chi_+(s, m_1, m_2)\chi_-(s, m_1, m_2))$
\rho(s, m_1, m_2) = \sqrt(\chi_+(s, m_1, m_2)\chi_-(s, m_1, m_2))