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


module mod_parameter_t_e

  implicit none

  integer(4), parameter :: IMAT_MAX = 4
  integer(4), parameter :: N_T_GRID_01 = 1000
  integer(4), parameter :: N_T_GRID_10 = 1000
  real(8), parameter :: E_EPS = 1.0e-5_8
  real(8), parameter :: CV_EPS = 1.0e-5_8
  
end module mod_parameter_t_e

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

module mod_variable_t_e

  implicit none

  real(8), allocatable :: t_grid(:,:)
  real(8), allocatable :: e_grid(:,:)
  
end module mod_variable_t_e


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

subroutine sub_create_t_e_table
  use mod_parameter_tillotson_eos
  use mod_parameter_t_e
  use mod_variable_t_e
  implicit none
  integer(4) :: mat,i
  real(8) :: debye,cv_lim
  real(8) :: fun_e_intg

  
  allocate(t_grid(IMAT_MAX,0:N_T_GRID_01+N_T_GRID_10))
  allocate(e_grid(IMAT_MAX,0:N_T_GRID_01+N_T_GRID_10))

  do mat=1,IMAT_MAX
    if (mat .eq. 1) then
       !---- Any Material ----
       debye=any_debye
       cv_lim=any_cv_lim
    elseif (mat .eq. 2) then
       !---- Murchison ----
       debye=murch_debye
       cv_lim=murch_cv_lim
    elseif (mat .eq. 3) then
       !---- Hard Ryugu ----
       debye=hryugu_debye
       cv_lim=hryugu_cv_lim
    elseif (mat .eq. 4) then
       !---- Soft Ryugu ----
       debye=sryugu_debye
       cv_lim=sryugu_cv_lim
    else
       write(*,*) '*************************************'
       write(*,*) '  MAT is invalid! Program stops!'
       write(*,*) '*************************************'
       stop
    endif

    t_grid(mat,0)=T_REF
    e_grid(mat,0)=0.0e0_8
    
    do i=1,N_T_GRID_01+N_T_GRID_10
       if(i > N_T_GRID_01) then
          t_grid(mat,i)=T_REF+real(N_T_GRID_01,8)+10.0e0_8*real(i-N_T_GRID_01,8)
       else
          t_grid(mat,i)=T_REF+real(i,8)
       endif       
       e_grid(mat,i)=fun_e_intg(debye,cv_lim,t_grid(mat,i),CV_EPS,E_EPS)
       !write(10,*) mat,i,t_grid(mat,i),e_grid(mat,i)
    enddo
    
 enddo
    
  return
end subroutine sub_create_t_e_table

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

real(8) function fun_e_intg(tdb,nk3,t_max,cv_eps,e_eps)
  use mod_parameter_tillotson_eos
  implicit none
  real(8), intent(IN) :: tdb,nk3,t_max,cv_eps,e_eps
  real(8) :: t_min
  integer(4) n,i
  real(8) :: dt,e,cv,t,e_old,e_err
  real(8) :: fun_cv_debye
  
  t_min=T_REF
  e_old=0.0e0_8

!  if(t_max <= t_min) then
!     fun_e_intg=0.0e0_8
!     go to 555
!  endif

if(t_max == t_min) then
fun_e_intg = 0.0e0_8
go to 555
endif
  
  n=100

999 continue
  
  dt=(t_max-t_min)/real(n,8)

  e=0.0e0_8
  
  do i=0,n
     t=t_min+dt*real(i,8)
     cv=fun_cv_debye(nk3,tdb,t,cv_eps)
     if(i == 0 .or. i == n) then
        e=e+cv*0.5e0_8*dt
     else
        e=e+cv*dt
     endif
  enddo

  e_err=abs((e-e_old)/e)

  !write(*,*) n,e_err,e
  
  if(e_err > e_eps) then
     n=n*2
     e_old=e
     go to 999
  endif
  !write(*,*) e/real(t_max-t_min),e

  fun_e_intg=e

555 continue
  
  return
end function fun_e_intg


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

real(8) function fun_cv_debye(nk3,tby,t,cv_eps)
  implicit none

  real(8), intent(IN) :: nk3,tby,t,cv_eps
  real(8), parameter :: NN = 100
  real(8) :: fun_debye
  real(8) :: a,b
  integer(4) :: n,i
  real(8) :: dx,intg,x
  real(8) :: intg_err
  real(8) :: intg_old
  real(8) :: pi

  if(t == 0.0e0_8) then
     fun_cv_debye = 0.0e0_8
     go to 555
  endif
  
  a=0.0e0_8
  b=tby/t

  if(b > 100.0e0_8) then
     pi=4.0e0_8*atan(1.0e0_8)
     fun_cv_debye=nk3*4.0e0_8*pi**4/5.0d0*(1.0e0_8/b)**3
     go to 555
  endif
  
  n=NN

  intg_old=0.0e0_8
  
999 continue
  dx=(b-a)/real(n,8)

  intg=0.0e0_8
  
  do i=1,n
     x=a+dx*real(i,8)
     if(i == n) then
        intg=intg+0.5e0_8*fun_debye(x)*dx
     else
        intg=intg+fun_debye(x)*dx
     endif
  enddo

  intg_err=abs((intg-intg_old)/intg)

  !write(*,*) n,intg_err
  
  if(intg_err > cv_eps) then
     n=n*2
     intg_old=intg
     go to 999
  endif
  
  fun_cv_debye=nk3*3.0e0_8*(t/tby)**3*intg

555 continue
  
  return
end function fun_cv_debye

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

real(8) function fun_debye(x)
  implicit none
  real(8), intent(IN) :: x

  fun_debye=x**4*exp(x)/(exp(x)-1.0e0_8)**2

  return

end function fun_debye