compute_post_shock Subroutine

   elemental subroutine compute_post_shock(eos, sgn, u0, p0, a0, ux, rx, px, ax, ss)
   !< Computing an unknown state `x` from a known state `0` when the two states are separated by a shock, given the velocity
   !< `ux`.
   !<
   !< The `sgn` dummy argument indicates if the shock propagates on `u-a (sgn=-1)` or `u+a (sgn=1)`.
   class(eos_object), intent(in)  :: eos        !< Equation of state.
   real(R8P),         intent(in)  :: sgn        !< Sign for distinguishing *left* (-1) from *right* (1) wave.
   real(R8P),         intent(in)  :: u0, p0, a0 !< Known state (speed, pressure and speed of sound).
   real(R8P),         intent(in)  :: ux         !< Unknown speed.
   real(R8P),         intent(out) :: rx, px, ax !< Unknown state (density, pressure and speed of sound).
   real(R8P),         intent(out) :: ss         !< Shock wave speed.
   real(R8P)                      :: M0         !< Relative Mach number of known state.
   real(R8P)                      :: x          !< Dummy variable.

   x   = 0.25_R8P * eos%gp1() * (ux - u0) / a0                     ! dummy variable
   M0  = x + sgn * sqrt(1.0_R8P + x * x)                           ! relative Mach number of known state
   x   = 1._R8P + 2._R8P * eos%g() * (M0 * M0 - 1._R8P) /eos%gp1() ! dummy variable (pressure ratio px/p0)

   ax = a0 * sqrt((eos%gp1() + eos%gm1() * x)/(eos%gp1() + eos%gm1() / x)) ! unknown speed of sound
   px = p0 * x                                                             ! unknown pressure
   rx = eos%density(pressure=px, speed_of_sound=ax)                        ! unknown density
   ss = u0 + a0 * M0                                                       ! shock wave speed
   endsubroutine compute_post_shock