Mathematica Formulation
Spherical parabolic coordinates are given by
\begin{equation}
x=\sigma\tau\cos\phi\\
\end{equation}
\begin{equation}
y=\sigma\tau\sin\phi\\
\end{equation}
\begin{equation}
z=\frac{1}{2}\left(\tau^{2}-\sigma^{2}\right)\\
\end{equation}
In these coordinates, the steady state, azimuthally symmetric hydrodynamics equations take the following form. The conservation of mass is
\begin{equation}
\frac{\partial}{\partial\sigma}\left(\sqrt{\sigma^{2}+\tau^{2}}\sigma\tau v_{\sigma}\right)+\frac{\partial}{\partial\tau}\left(\sqrt{\tau^{2}+\sigma^{2}}\sigma\tau v_{\tau}\right)=0\\
\end{equation}
The conservation of entropy is
\begin{equation}
v_{\sigma}\frac{\partial s}{\partial\sigma}+v_{\tau}\frac{\partial s}{\partial\tau}=0\\
\end{equation}
where \(s=\ln p-\gamma\ln\rho\) is the specific entropy.
The conservation of momentum in each direction is
\begin{equation}
v_{\sigma}\frac{\partial v_{\sigma}}{\partial\sigma}+v_{\tau}\frac{\partial v_{\sigma}}{\partial\tau}+\frac{v_{\tau}\left(v_{\sigma}\tau-v_{\tau}\sigma\right)}{\tau^{2}+\sigma^{2}}+\frac{1}{\rho}\frac{\partial p}{\partial\sigma}=0\\
\end{equation}
\begin{equation}
v_{\sigma}\frac{\partial v_{\tau}}{\partial\sigma}+v_{\tau}\frac{\partial v_{\tau}}{\partial\tau}-\frac{v_{\sigma}\left(v_{\sigma}\tau-v_{\tau}\sigma\right)}{\tau^{2}+\sigma^{2}}+\frac{1}{\rho}\frac{\partial p}{\partial\tau}=0\\
\end{equation}
One of the momentum equations can be replaced by Bernoulli’s equation
\begin{equation}
\frac{1}{2}v_{in}^{2}=\frac{1}{2}v_{\sigma}^{2}+\frac{1}{2}v_{\tau}^{2}+\frac{1}{\gamma-1}c^{2}\\
\end{equation}
where \(c\) is the speed of sound. We found it most convenient to express the equations in terms of the density \(\rho\), the sound speed \(c\), and the two Mach numbers \(m_{\sigma,\tau}=\frac{v_{\sigma,\tau}}{c}\)
We assume that the shock front coincides with a curve on which \(\sigma=\sigma_{s}\) is constant. The Rankine Hugoniot boundary conditions at the shock fronts are
\begin{equation}
\rho_{s}=\frac{\gamma+1}{\gamma-1}\rho_{a}\\
\end{equation}
\begin{equation}
c_{s}=\frac{\sqrt{2\gamma\left(\gamma-1\right)}}{\gamma+1}v_{in}\frac{\sigma_{s}}{\sqrt{\tau^{2}+\sigma_{s}^{2}}}\\
\end{equation}
\begin{equation}
m_{\sigma s}=-\sqrt{\frac{\gamma-1}{2\gamma}}\\
\end{equation}
\begin{equation}
m_{\tau s}=\frac{\gamma+1}{\sqrt{2\gamma\left(\gamma-1\right)}}\frac{\tau}{\sigma_{s}}\\
\end{equation}
So we have a well defined boundary value problem. However, since it involves two coordinates, our approach so far has no benefit over the hydrodynamic equations in cylindrical coordinates. At very high values of \(\tau\gg\sigma_{s}\), the equations can be simplified. We assume that the variables vary with \(\tau\) as the shocked values do, i.e. \(\rho\) is independent of \(\tau\), \(c\propto\frac{1}{\tau}\), \(m_{\sigma}\) is independent of \(\tau\) and \(m_{\tau}\propto\tau\). Using this approximation, \(\tau\) can be eliminated from the hydrodynamic equations, and the problem reduces to a set of ordinary differential equation in \(\sigma\). These equation can be numerically integrated, and we can obtain curves for \(\rho\), \(\tau\cdot c\), \(m_{\sigma}\) and \(\frac{m_{\tau}}{\tau}\). We performed the complete analysis in a mathematica notebook.
