hugoniot_star_state Subroutine

subroutine hugoniot_star_state(sgnBC, Pstar, Ro_LR, p_LR, u_LR, c_LR, st)
  real(dp),         intent(in)  :: sgnBC, Pstar, Ro_LR, p_LR, u_LR, c_LR
  type(hug_star_t), intent(out) :: st

  integer,  parameter :: NbreMaxIte = 200
  real(dp), parameter :: MachError = 1.0e-15_dp, Tol_brent = 1.0e-10_dp

  real(dp) :: DeltaP, v_LR
  real(dp) :: e_LR, ev_LR, ep_LR
  real(dp) :: v, v_pred, v_lo, v_hi, Psi
  real(dp) :: e_v, ev_v, ep_v
  real(dp) :: Funct, DerFunct, D, tol_H, droeint_dro
  integer  :: ite
  logical  :: have_lo, have_hi
  type(hug_brent_t) :: hug_obj

  st%ok     = .false.
  st%n_iter = 0

  DeltaP = Pstar - p_LR
  v_LR   = 1.0_dp / Ro_LR

  droeint_dro = droeint_droP(Ro_LR, p_LR, e=e_LR, dedv=ev_LR, dedp=ep_LR)

  v_pred = v_LR - DeltaP * v_LR**2 / c_LR**2
  if (v_pred <= 0.0_dp) v_pred = v_LR * 0.5_dp
  v = v_pred

  hug_obj%Pstar = Pstar
  hug_obj%p_LR  = p_LR
  hug_obj%v_LR  = v_LR
  hug_obj%e_LR  = e_LR

  v_lo = 0.0_dp; v_hi = 0.0_dp
  have_lo = .false.; have_hi = .false.

  tol_H = 1.0e-8_dp * (Pstar + p_LR) * v_LR

  do ite = 1, NbreMaxIte
    droeint_dro = droeint_droP(1.0_dp/v, Pstar, e=e_v, dedv=ev_v, dedp=ep_v)
    Funct   = (Pstar + p_LR)*(v - v_LR) + 2.0_dp*(e_v - e_LR)
    DerFunct = (Pstar + p_LR) + 2.0_dp*ev_v
    if (Funct < 0.0_dp) then
      if (.not. have_lo .or. v > v_lo) then; v_lo = v; have_lo = .true.; endif
    else
      if (.not. have_hi .or. v < v_hi) then; v_hi = v; have_hi = .true.; endif
    endif
    if (abs(Funct) > tol_H) then
      v = v - Funct/DerFunct
      if (v <= 0.0_dp) exit
    else
      st%n_iter = ite; exit
    endif
    if (ite == NbreMaxIte) then
      if (have_lo .and. have_hi) then
        v = zero(hug_obj, v_lo, v_hi, MachError, Tol_brent)
        droeint_dro = droeint_droP(1.0_dp/v, Pstar, e=e_v, dedv=ev_v, dedp=ep_v)
      else
        call set_error('hugoniot_star_state: Newton did not converge, no bracket')
        return
      endif
      exit
    endif
  enddo

  if (v <= 0.0_dp) then
    if (.not. (have_lo .and. have_hi)) then
      call set_error('hugoniot_star_state: cannot bracket root')
      return
    endif
    v = zero(hug_obj, v_lo, v_hi, MachError, Tol_brent)
    droeint_dro = droeint_droP(1.0_dp/v, Pstar, e=e_v, dedv=ev_v, dedp=ep_v)
  endif

  Psi = (v_LR - v) / DeltaP
  if (Psi <= 0.0_dp) then
    call set_error('hugoniot_star_state: Psi = (v_LR - v*)/DeltaP <= 0')
    return
  endif

  st%vstar = v
  st%K     = sqrt(Psi)
  st%W     = 1.0_dp / st%K
  st%ustar = u_LR + sgnBC * DeltaP * st%K
  st%ok    = .true.

  D = (Pstar + p_LR) + 2.0_dp * ev_v
  st%dvstar_dPstar = -((v - v_LR) + 2.0_dp * ep_v) / D
  st%dustar_dPstar = sgnBC * (st%K + (-st%dvstar_dPstar - Psi) / (2.0_dp * st%K))
  st%dvstar_dv_LR  = ((Pstar + p_LR) + 2.0_dp * ev_LR) / D
  st%dvstar_dp_LR  = -((v - v_LR) - 2.0_dp * ep_LR) / D
  st%ev_v          = ev_v

end subroutine hugoniot_star_state