add debug checks for cons2prim

src/projects/GRHD/cons2prim_implementations.jl

 @with_signature function cons2prim(equation::Equation{<:HotEquationOfState})
   @accepts D, Sr, tau, p, r, B
-  S2   = (Sr/B)^2
-  @unpack eos, cons2prim_newtonsolver, cons2prim_bisection = equation
-  p, converged = find_pressure_insane(D, S2, tau, p, r, eos,
-                                      cons2prim_newtonsolver, cons2prim_bisection)
-  if !converged && eltype(equation) == HybridEoS
-    rho, success = find_density_insane(D, S2, p, r, equation, eos,
-                                       cons2prim_newtonsolver, cons2prim_bisection)
-    if !success
-      error()
-    end
-    p    = eos(Pressure, Density, rho)
-    vr   = Formulae.vr(B, D, Sr, tau, p)
-    eps  = Formulae.ϵ(B, D, Sr, tau, p)
-  else
-    rho  = Formulae.ρ(B, D, Sr, tau, p)
-    vr   = Formulae.vr(B, D, Sr, tau, p)
-    eps  = Formulae.ϵ(B, D, Sr, tau, p)
-  end
+  # S2   = (Sr/B)^2
+  # @unpack eos, cons2prim_newtonsolver, cons2prim_bisection = equation
+  # p, converged = find_pressure_insane(D, S2, tau, p, r, eos,
+  #                                     cons2prim_newtonsolver, cons2prim_bisection)
+  # if !converged && eltype(equation) == HybridEoS
+  #   rho, success = find_density_insane(D, S2, p, r, equation, eos,
+  #                                      cons2prim_newtonsolver, cons2prim_bisection)
+  #   if !success
+  #     error()
+  #   end
+  #   p    = eos(Pressure, Density, rho)
+  #   vr   = Formulae.vr(B, D, Sr, tau, p)
+  #   eps  = Formulae.ϵ(B, D, Sr, tau, p)
+  # else
+  #   rho  = Formulae.ρ(B, D, Sr, tau, p)
+  #   vr   = Formulae.vr(B, D, Sr, tau, p)
+  #   eps  = Formulae.ϵ(B, D, Sr, tau, p)
+  # end
+  rho, vr, eps, p = cons2prim_kastaun((D, Sr, tau, p, r, B), equation)
   @returns rho, vr, eps, p
 end
 
@@ -134,10 +135,11 @@ end
 @with_signature function cons2prim_kastaun(equation::AbstractEquation)
   @accepts D, Sr, tau, p, r, B
 
-  @unpack eos = equation
-  rhomin, rhomax = 1e-10, Inf
-  fn_epsmin, fn_epsmax = rho -> -(1-1e-6), rho -> Inf
+  @unpack eos, atmosphere = equation
+  rhomin, rhomax = atmosphere.atm_density * atmosphere.threshold, Inf
+  fn_epsmin, fn_epsmax = rho -> -(1-1e-10), rho -> Inf
   h0 = 1e-3
+  rhoatm = atmosphere.atm_density
 
   q = tau / D
   S2 = (Sr/B)^2
@@ -150,86 +152,127 @@ end
   # What is not a question but rather W^hat, cf. (40)
   function fn_What(mu)
     _vhat = min(mu * s, v0)
+    if _vhat < 0
+      println((; q, _vhat, s, mu, v0))
+      error()
+    end
     return 1 / sqrt(1 - _vhat^2)
   end
 
   # setup bracket
   mu_a = 0.0
-  # with no EM fields we can find the root of (49) analytically
+  # with no EM fields we can find the root of (49) analytically, because s = const.
   mu_plus = 1 / sqrt(h02 + s2)
   mu_b = r < h0 ? 1 / h0 : mu_plus
 
   What_b = fn_What(mu_b)
-  if D / What_b > rhomax || D < rhomin
-    error("Invalid state encountered!")
+  if D / What_b > rhomax
+    error("Invalid state encountered at @ r = $(r)!")
   end
+  if D < rhomin
+    # @info "Invalid state encountered at @ r = $(r)!"
+    # @info "Atmosphering ..."
+    # println((; D, rhomax, rhomin))
+    rho = rhoatm
+    eps = eos.cold_eos(InternalEnergy, Density, rho)
+    vr = 0.0
+    p = eos.cold_eos(Pressure, Density, rho)
+    # p = eos(Pressure, Density, rho, InternalEnergy, eps)
+  else
 
-  # do we have to adjust the bracket?
-  if D / What_b < rhomax < D
-    pblm_rhomax = Roots.ZeroProblem(mu -> D / fn_What(mu) - rhomax, (mu_a, mu_b))
-    mu_b = Roots.solve(pblm_rhomax, Roots.A42())
-    if isnan(mu_b)
-      error("Failed to correct upper interval bound")
+    # do we have to adjust the bracket?
+    debug = false
+    # if D / What_b < rhomax < D
+    if rhomax < D
+      pblm_rhomax = Roots.ZeroProblem(mu -> D / fn_What(mu) - rhomax, (mu_a, mu_b))
+      mu_b = Roots.solve(pblm_rhomax, Roots.A42())
+      if isnan(mu_b)
+        error("Failed to correct upper interval bound")
+      end
+      debug = true
     end
-  end
-  if D / What_b < rhomin < D
-    pblm_rhomin = Roots.ZeroProblem(mu -> D / fn_What(mu) - rhomin, (mu_a, mu_b))
-    mu_a = Roots.solve(pblm_rhomin, Roots.A42())
-    if isnan(mu_a)
-      error("Failed to correct lower interval bound")
+    # if D / What_b < rhomin < D
+    if D / What_b < rhomin
+      mu_rng = collect(LinRange(mu_a, mu_b, 1000))
+      fn = @. D / fn_What(mu_rng) - rhomin
+      data = hcat(mu_rng, fn)
+      # display(data)
+      # display(D / What_b)
+      # display(D)
+      # display(rhomin)
+      pblm_rhomin = Roots.ZeroProblem(mu -> D / fn_What(mu) - rhomin, (mu_a, mu_b))
+      mu_b = Roots.solve(pblm_rhomin, Roots.A42())
+      if isnan(mu_b)
+        error("Failed to correct lower interval bound")
+      end
+      println("New interval: ($mu_a, $mu_b) @ r = $r")
+      debug = false
+    end
+
+    function master_fn(mu)
+      vhat = min(mu * s, v0)
+      What = 1 / sqrt(1 - vhat^2)
+      rhohat = D / What
+      eps_plus_one_over_What = What * (q - mu * s2)
+      epshat = eps_plus_one_over_What + vhat^2 * What^2 / (1 + What)
+      # enforce EOS bounds
+      _epsmin, _epsmax = fn_epsmin(rhohat), fn_epsmax(rhohat)
+      if rhohat < rhomin || epshat < _epsmin
+        rhohat, epshat, What = rhomin, _epsmin
+      elseif rhohat > rhomax || epshat > _epsmax
+        rhohat, epshat = rhomax, _epsmax
+      end
+      # compute thermodynamic vars
+      phat = eos(Pressure, Density, rhohat, InternalEnergy, epshat)
+      ahat = phat / (rhohat * (1 + epshat))
+      hhat = eos(Enthalpy, Density, rhohat, InternalEnergy, epshat)
+      # compute master function
+      # extracted from (35), see hint in paragraph after (46)-(48)
+      nu_A = (1 + ahat) * eps_plus_one_over_What
+      nu_B = (1 + ahat) * (1 + q - mu * s2)
+      nu = max(nu_A, nu_B)
+      f = mu - 1 / (nu + s2 * mu)
+      return f
     end
-  end
 
+    if debug
+      mu_rng = collect(LinRange(mu_a, mu_b, 1000))
+      fn = @. master_fn(mu_rng)
+      data = hcat(mu_rng, fn)
+      display(data)
+    end
 
-  function master_fn(mu)
-    vhat = min(mu * s, v0)
-    What = 1 / sqrt(1 - vhat^2)
-    rhohat = D / What
-    eps_plus_one_over_What = What * (q - mu * s2)
-    epshat = eps_plus_one_over_What + vhat^2 * What^2 / (1 + What)
+    mf_a = master_fn(mu_a)
+    mf_b = master_fn(mu_b)
+    if mf_a * mf_b >= 0
+      println((; D, Sr, tau, p, r, B, mf_a, mf_b, mu_a, mu_b))
+      error("Master function is ill conditioned, @ r = $(r)!")
+    end
+    pblm = Roots.ZeroProblem(master_fn, (mu_a, mu_b))
+    mu = Roots.solve(pblm, Roots.A42() #= =^= algorithm 748 =#)
+    if isnan(mu)
+      error("failed to find root")
+    end
+
+    # evaluate primitives
+    v = min(mu * s, v0)
+    W = 1 / sqrt(1 - v^2)
+    rho = D / W
+    eps_plus_one_over_W = W * (q - mu * s2)
+    eps = eps_plus_one_over_W + v^2 * W^2 / (1 + W)
     # enforce EOS bounds
-    _epsmin, _epsmax = fn_epsmin(rhohat), fn_epsmax(rhohat)
-    if rhohat < rhomin || epshat < _epsmin
-      rhohat, epshat = rhomin, _epsmin
-    elseif rhohat > rhomax || epshat > _epsmax
-      rhohat, epshat = rhomax, _epsmax
+    epsmin, epsmax = fn_epsmin(rho), fn_epsmax(rho)
+    if rho < rhomin || eps < epsmin
+      rho, eps = rhomin, epsmin
+    elseif rho > rhomax || eps > epsmax
+      rho, eps = rhomax, epsmax
     end
     # compute thermodynamic vars
-    phat = eos(Pressure, Density, rhohat, InternalEnergy, epshat)
-    ahat = phat / (rhohat * (1 + epshat))
-    hhat = eos(Enthalpy, Density, rhohat, InternalEnergy, epshat)
-    # compute master function
-    # extracted from (35), see hint in paragraph after (46)-(48)
-    nu_A = (1 + ahat) * eps_plus_one_over_What
-    nu_B = (1 + ahat) * (1 + q - mu * s2)
-    nu = max(nu_A, nu_B)
-    f = mu - 1 / (nu + s2 * mu)
-    return f
-  end
+    p = eos(Pressure, Density, rho, InternalEnergy, eps)
 
-  pblm = Roots.ZeroProblem(master_fn, (mu_a, mu_b))
-  mu = Roots.solve(pblm, Roots.A42() #= =^= algorithm 748 =#)
-  if isnan(mu)
-    error("failed to find root")
-  end
+    vr = Sr*mu/D
 
-  # evaluate primitives
-  v = min(mu * s, v0)
-  W = 1 / sqrt(1 - v^2)
-  rho = D / W
-  eps_plus_one_over_W = W * (q - mu * s2)
-  eps = eps_plus_one_over_W + v^2 * W^2 / (1 + W)
-  # enforce EOS bounds
-  epsmin, epsmax = fn_epsmin(rho), fn_epsmax(rho)
-  if rho < rhomin || eps < epsmin
-    rho, eps = rhomin, epsmin
-  elseif rho > rhomax || eps > epsmax
-    rho, eps = rhomax, epsmax
   end
-  # compute thermodynamic vars
-  p = eos(Pressure, Density, rho, InternalEnergy, eps)
-
-  vr = Sr*mu/D
 
   @returns rho, vr, eps, p
 end