subroutine compute_primitive_from_conserved(U, prim, error)

  ! This function computes primitive variables from conserved variables, using Newton iterations.
  ! A constant-Gamma equation of state is employed.
  ! The quantity "W" is the Lorentz factor (W letter is used to avoid confusion with the adiabatic 
  ! gas constant "Gamma")
  ! The rest of the nomenclature follows the equations in the paper by Mignone and Bodo, 
  ! "An HLLC Riemann solver for relativistic flows – I. Hydrodynamics", 2005.

  implicit none

  real(kind=8), dimension(4), intent(in)  :: U
  real(kind=8), dimension(4), intent(out) :: prim
  real(kind=8), intent(out) :: error

  real(kind=8) :: D, mx, my, E, m2
  real(kind=8) :: W, rhoh, rho, ux, uy, P
  real(kind=8) :: c, Gam

  real(kind=8) :: H, dHdP
  real(kind=8) :: Pnew, tol
  integer      :: itermax, iter_count

  c   = 1.0d0       ! Use units where c = 1
  Gam = 4.0d0/3.0d0 ! Constant Gamma = 4/3

  ! Extract conserved variables
  D  = U(1)
  mx = U(2)
  my = U(3)
  E  = U(4)

  m2 = mx*mx + my*my ! Working variable

  ! Initialize Newton solver
  tol        = 1.0d-15   ! Accepted error between iterations
  itermax    = 1000      ! Maxium number of iterations
  iter_count = 0         ! Initialize

  P      = 1.0d0 ! Init
  error  = 1.0d0 ! init

  ! Iterate until convergence
  do while (error .gt. tol)
  
    print*, 'Current P: ', P, '  error: ', error
  
    ! Function to be minimized...
    H = D*c*c*sqrt(P*P + 2.0d0*E*P + E*E - m2/c/c) + 1.0d0/(Gam-1.0d0)*P*P &
       + E*(2.0d0-Gam)/(Gam-1.0d0)*P + m2/c/c - E*E
    
    ! ... and its derivative
    dHdP = 2.0d0/(Gam-1.0d0)*P + E*(2.0d0-Gam)/(Gam-1.0d0) &
         + D*c*c*(P+E)/sqrt(P*P + 2.0d0*E*P + E*E - m2/c/c)
  
    Pnew  = P - H/dHdP    ! New pressure
    error = abs(Pnew - P) ! Compute error
    P     = Pnew          ! Update P
  
    iter_count = iter_count + 1
    if (iter_count .ge. itermax) then
      print*, "ATTENTION!"
      print*, "Iteration limit exceeded while computing primitive variables!"
      exit ! Exit do loop
    end if
  
  end do

  ! Compute other quantities
  W    = 1.0d0/sqrt(1.0d0 - m2/(c*c*(E+P)*(E+P))) ! Lorentz factor
  rho  = D/W
  rhoh = rho*c*c + P*Gam/(Gam-1.0d0)
  ux   = mx/(rhoh*W*W)
  uy   = my/(rhoh*W*W)

  ! Compose array of primitive variables
  prim(1) = rho
  prim(2) = ux
  prim(3) = uy
  prim(4) = P

end subroutine

! ----------------------------------------
! MAIN PROGRAM
! ----------------------------------------

program test_inverse

! This program tests the computation of primitive variables rho, ux, uy, P
! from conserved variables D, mx, my, E.
! See the subroutine above for more details.
! 
! The idea here is that:
! 1) we start by some known primitive variables, that we set arbitrarily
! 2) we compute the corresponding conserved variables D, mx, my, E (this can
!    be done easily by hand, of course)
! 3) we invert these conserved variables numerically, using the Newton method
!    implemented in the subroutine

implicit none

real(kind=8) :: c, Gam, g 
real(kind=8) :: rho, ux, uy, P ! Primitive variables

real(kind=8) :: D, rhoh, mx, my, E ! Conserved variables
real(kind=8) :: u2, m2, error ! Working variables

real(kind=8), dimension(4) :: U, prim

integer :: i

c   = 1.0d0
Gam = 4.0d0/3.0d0

! Set primitive variables
rho = 2.0d0
ux  = 0.3d0*c
uy  = 0.0d0
P   = 8.843d0

! Create conserved variables
u2 = ux*ux + uy*uy
g  = 1.0d0/sqrt(1.0d0 - u2/(c*c))

D    = g*rho
rhoh = rho*c*c + P*Gam/(Gam - 1.0d0)
mx   = rhoh*g*g*ux
my   = rhoh*g*g*uy
E    = rhoh*g*g - P

print*, "TARGET PRIMITIVE VARIABLES: "
print*, "   rho: ", rho
print*, "   ux:  ", ux
print*, "   uy:  ", uy
print*, "   P:   ", P
print*, " "
print*, " gamma: ", g
print*, " "
print*, "CONSERVED VARIABLES:"
print*, "   D:  ", D
print*, "   mx: ", mx
print*, "   my: ", my
print*, "   E:  ", E

! Now, we have the conserved variables, we want to recompute the
! primitive variables using Newton iterations

U(1) = D
U(2) = mx
U(3) = my
U(4) = E

call compute_primitive_from_conserved(U, prim, error)

rho = prim(1)
ux  = prim(2)
uy  = prim(3)
P   = prim(4)

print*, 'DBDBDB', prim

print*, " "
print*, "NUMERICAL PRIMITIVE VARIABLES:"
print*, "   rho: ", rho
print*, "   ux:  ", ux
print*, "   uy:  ", uy
print*, "   P:   ", P, "  (final error: ", error, ")"

end program