!program main_test
!  implicit none
!
!  call sub_create_rho_ec_table
!
!  stop
!end program main_test

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------


module mod_parameter_ec

  implicit none

  integer(4), parameter :: IMAT_MAX = 4
  integer(4), parameter :: N_RHO_GRID = 1000
  real(8),    parameter :: ALP_RHO = 1.003d0
  real(8),    parameter :: EPS_EC = 1.0d-5
  
end module mod_parameter_ec


!-----------------------------------------------------------------------
!-----------------------------------------------------------------------


module mod_variable_ec_tillotson

  implicit none

  real(8) :: rho0,a,b,ap,bp,e0
  
end module mod_variable_ec_tillotson


!-----------------------------------------------------------------------
!-----------------------------------------------------------------------


module mod_variable_ec

  implicit none

  real(8), allocatable :: rho_grid(:,:)
  real(8), allocatable :: ec_grid(:,:)
  
end module mod_variable_ec


!-----------------------------------------------------------------------
!-----------------------------------------------------------------------


subroutine sub_create_rho_ec_table
  use mod_parameter_tillotson_eos
  use mod_parameter_ec
  use mod_variable_ec
  use mod_variable_ec_tillotson
  implicit none

  real(8) :: rho,ec
  integer(4) :: mat,i

  allocate(rho_grid(IMAT_MAX,0:N_RHO_GRID))
  allocate(ec_grid(IMAT_MAX,0:N_RHO_GRID))
  
  do mat=1,IMAT_MAX
    if (mat .eq. 1) then
       !---- Any Material ----
       rho0=any_rho0
       a=any_a
       b=any_b
       ap=any_ap
       bp=any_bp
       e0=any_e0
    elseif (mat .eq. 2) then
       !---- Murchison ----
       rho0=murch_rho0
       a=murch_a
       b=murch_b
       ap=murch_ap
       bp=murch_bp
       e0=murch_e0
    elseif (mat .eq. 3) then
       !---- Hard Ryugu ----
       rho0=hryugu_rho0
       a=hryugu_a
       b=hryugu_b
       ap=hryugu_ap
       bp=hryugu_bp
       e0=hryugu_e0
    elseif (mat .eq. 4) then
       !---- Soft Ryugu ----
       rho0=sryugu_rho0
       a=sryugu_a
       b=sryugu_b
       ap=sryugu_ap
       bp=sryugu_bp
       e0=sryugu_e0
    else
       write(*,*) '*************************************'
       write(*,*) '  MAT is invalid! Program stops!'
       write(*,*) '*************************************'
       stop
    endif

    rho=rho0

    do i=0,N_RHO_GRID
       call sub_calc_ec(rho,ec)
       rho_grid(mat,i)=rho
       ec_grid(mat,i)=ec
       !if (mat == 2) write(10,*) rho,ec
       rho=rho*ALP_RHO
    enddo

!    write(10,*) ' '

 enddo

 return
end subroutine sub_create_rho_ec_table


!=======================================================================
!=======================================================================


subroutine sub_calc_ec(rho,ec)
  use mod_parameter_ec
  use mod_variable_ec_tillotson
  implicit none

  real(8), intent(IN)  :: rho
  real(8), intent(OUT) :: ec

  integer(4) :: nn,n
  real(8) :: drho,ec_new,ec_old,rho_tmp,decdrho
  real(8) :: ec_err
  real(8) :: p_till

  if(rho <= rho0) then
     ec=0.0d0
     goto 555
  endif

  nn=100

  ec_old=1.0e30_8

  do
     drho=(rho-rho0)/real(nn,8)
     ec_new=0.0e0_8
     do n=0,nn-1
        rho_tmp=rho0+drho*real(n,8)
        decdrho=p_till(rho_tmp,ec_new)/rho_tmp/rho_tmp
        ec_new = ec_new + decdrho * drho
     enddo
     ec_err=abs((ec_old-ec_new)/ec_new)
     !write(30,*) rho,nn,ec_err
     if(ec_err < EPS_EC) exit
     nn=2*nn
     ec_old=ec_new
  enddo

  ec=ec_new

  555 continue

  return
end subroutine sub_calc_ec


!=======================================================================
!=======================================================================


real(8) function p_till(rho,e)

  use mod_variable_ec_tillotson

  implicit none

!---- Arguments in This Subroutine -----------------------------
  real(8), intent(IN)  ::  rho       ! Density [kg/m^3]
  real(8), intent(IN)  ::  e         ! Internal Energy [J/kg]

!---- Variables Used in This Subroutine ------------------------
  real(8) :: p
  real(8) :: eta
  real(8) :: um
!---------------------------------------------------------------

  eta=rho/rho0
  um=eta-1.0d0

  p=(a+b/(e/e0/eta/eta+1.0d0))*rho*e
  p=p+ap*um+bp*um*um
  if (p .lt. 0.0d0) p=0.0d0

  p_till=p
      
  return
end function p_till