!! This simple procedure is necessary to generate the h5 table



program generate_water_table

  use aneos_access

  use ideal_water

  implicit none

  

  integer :: i, j, NRho, NT, outid, outid2, ios

  double precision :: logRho, logT, lnRho, lnT, Rho, T, Fvect(num_free_derivs)

  double precision, pointer, dimension(:) :: logRhovect, logTvect

  double precision, pointer, dimension(:,:) :: Fr, dfdrho, dfdt, d2fdrho2, &

       d2fdt2, d2fdrhodt, d3fdrho2dt, d3fdrhodt2, d4fdrho2dt2

  double precision :: logRhoLow, logRhoHi, dlogRho

  double precision :: logTLow, logTHi, dlogT

  character (len=256) :: filename

  double precision :: f_val, df_drho, df_dt, d2f_drho2, d2f_dt2, d2f_drhodt, &

       d3f_drho3, d3f_drho2dt, d3f_drhodt2, d3f_dt3, d4f_drho2dt2



  NRho = 81

  NT = 31

  

  allocate(Fr(NRho, NT))

  allocate(dfdrho(NRho, NT))

  allocate(dfdt(NRho, NT))

  allocate(d2fdrho2(NRho, NT))

  allocate(d2fdt2(NRho, NT))

  allocate(d2fdrhodt(NRho, NT))

  allocate(d3fdrho2dt(NRho, NT))

  allocate(d3fdrhodt2(NRho, NT))

  allocate(d4fdrho2dt2(NRho, NT))



  allocate(logRhovect(NRho))

  allocate(logTvect(NT))



  logRhoLow = -5.

  logRhoHi = 3.

  logTLow = 2.8

  logTHi = 6.8

  

  dlogRho = (logRhoHi - logRhoLow) / (NRho - 1)

  dlogT = (logTHi - logTLow) / (NT - 1)

  

  !! Write out both formatted and unformatted data

  !! I am trying to shift to the unformatted approach since it

  !!  1. takes up less space

  !!  2. Doesn't introduce more numerical errors 

  outid = 51

  outid2 = 52

  open(unit=outid,file='data/wateridealfe.bin', form='unformatted')  

  open(unit=outid2,file='data/wateridealfe.dat')

  

  do i=1,NRho

     do j=1,NT

        logRho = logRhoLow + (i-1) * dlogRho

        logT = logTLow + (j-1) * dlogT

        logRhovect(i) = logRho

        logTvect(j) = logT

        

        Rho = 10.d0**logRho

        T = 10.d0**logT

        

        lnRho = log(Rho)

        lnT = log(T)

        

        call get_IG_EOS(lnRho, lnT, f_val, df_drho, df_dt, &

             d2f_drho2, d2f_dt2, d2f_drhodt, &

             d3f_drho3, d3f_drho2dt, d3f_drhodt2, d3f_dt3, &

             d4f_drho2dt2)



        Fvect(a_val) = f_val

        Fvect(a_dRho) = df_drho

        Fvect(a_dT) = df_dt

        Fvect(a_dRho2) = d2f_drho2

        Fvect(a_dT2) = d2f_dt2

        Fvect(a_dRhodT) = d2f_drhodt

        Fvect(a_dRho2dT) = d3f_drho2dt

        Fvect(a_dRhodT2) = d3f_drhodt2

        Fvect(a_dRho2dT2) = d4f_drho2dt2

        

        write(outid) lnRho, lnT, Fvect

        write(outid2,*) lnRho, lnT, Fvect



        Fr(i,j) = Fvect(a_val)

        dfdrho(i,j) = Fvect(a_dRho)

        dfdt(i,j) = Fvect(a_dT)

        d2fdrho2(i,j) = Fvect(a_dRho2)

        d2fdt2(i,j) = Fvect(a_dT2)

        d2fdrhodt(i,j) = Fvect(a_dRhodT)

        d3fdrho2dt(i,j) = Fvect(a_dRho2dT)

        d3fdrhodt2(i,j) = Fvect(a_dRhodT2)

        d4fdrho2dt2(i,j) = Fvect(a_dRho2dT2)

        

     enddo

  enddo



  filename = "h5tdata/fr_logRhovect.data"

  open(unit=19,file=trim(filename),status='replace', &

         iostat=ios,action='write',form='formatted')     

  write(19,'(e20.6)') logRhovect

  close(unit=19)

  

  filename = "h5tdata/fr_logTvect.data"

  open(unit=19,file=trim(filename),status='replace', &

         iostat=ios,action='write',form='formatted')     

  write(19,'(e20.6)') logTvect

  close(unit=19)



  filename = "h5tdata/fr.data"

  call write_matrix(Fr, filename)

  

  filename = "h5tdata/dfdrho.data"

  call write_matrix(dfdrho, filename)

  

  filename = "h5tdata/dfdt.data"

  call write_matrix(dfdt, filename)  



  filename = "h5tdata/d2fdrho2.data"

  call write_matrix(d2fdrho2, filename)

  

  filename = "h5tdata/d2fdt2.data"

  call write_matrix(d2fdt2, filename)

  

  filename = "h5tdata/d2fdRhodT.data"

  call write_matrix(d2fdrhodt, filename)

  

  filename = "h5tdata/d3fdRho2dT.data"

  call write_matrix(d3fdrho2dt, filename)

  

  filename = "h5tdata/d3fdRhodT2.data"

  call write_matrix(d3fdrhodt2, filename)

  

  filename = "h5tdata/d4fdRho2dT2.data"

  call write_matrix(d4fdrho2dt2, filename)

  

  close(unit=outid)

  close(unit=outid2)

  deallocate(Fr)

  deallocate(dfdrho)

  deallocate(dfdt)

  deallocate(d2fdrho2)

  deallocate(d2fdt2)

  deallocate(d2fdrhodt)

  deallocate(d3fdrho2dt)

  deallocate(d3fdrhodt2)

  deallocate(d4fdrho2dt2)

  deallocate(logRhovect)

  deallocate(logTvect)



  





contains

  subroutine write_matrix(matrix, filename)

    double precision, intent(in), dimension(:,:), pointer :: matrix

    character (len=256) :: filename

    

    open(unit=19,file=trim(filename),status='replace', &

         iostat=ios,action='write',form='formatted')

    write(19,'(e20.6)') matrix

    close(unit=19)

  end subroutine write_matrix



end program generate_water_table