Solution (source code)

= Solution

The radial equation can be written
$$
\frac{d}{dr}\left(\frac{u^2}{2}
+c_s^2\ln\rho-\frac{GM}{r}\right)=0.
$$
Thus the <Bernoulli function> is constant on each <streamline>:
$$
\boxed{\mathcal B=\frac{u^2}{2}+c_s^2\ln\rho+\Phi},
\qquad
\Phi=-\frac{GM}{r}.
$$
The boundary conditions $u\to0$ and $\rho\to\rho_0$ as $r\to\infty$ fix $\mathcal B=c_s^2\ln\rho_0$.