! See copyright notice in the COPYRIGHT file.
!> This module contains an interface for external C++ code to compute
!! liquid mixture property like thermodynamic factor and 
!! Maxwell-Stefan Diffusivity coefficients
module mus_eNRTL_module
  use, intrinsic :: iso_c_binding, only: c_ptr, c_int, c_double, c_float,      &
    &                                    c_int_least8_t, c_int_least64_t,      &
    &                                    c_char, c_loc, c_bool

  use env_module,                  only: rk  
  
  implicit none

  private

  public :: mus_init_eNRTL
  public :: mus_calc_thermFactor
  public :: mus_calc_MS_DiffMatrix

  !> This function initialize eNRTL model by loading liquid mixture
  !! property from filename
  interface
    function init_enrtl_loc( filename, nSpc ) bind(c, name='init_enrtl')
    !function init_enrtl_f90( filename ) bind(c, name='_Z10init_enrtlPKc')
      use, intrinsic :: iso_c_binding
      character(kind=c_char), dimension(*) :: filename
      integer(kind=c_int), intent(out) :: nSpc
      !> Result, indicating the status of encode
      integer(kind=c_int) :: init_enrtl
    end function init_enrtl_loc
  end interface  

  !> This routine calculates thermodynamic factor for given mole_frac
  !! of all species
  interface
    subroutine calc_therm_factor_loc( nSpc, Temp, Press, Mole_frac,            &
      & Therm_factors ) bind(c, name='calc_therm_factor_C')
      use, intrinsic :: iso_c_binding
      integer(kind=c_int), value, intent(in) :: nSpc
      real(kind=c_double), value, intent(in) :: Temp
      real(kind=c_double), value, intent(in) :: Press
      real(kind=c_double), dimension(*), intent(in) :: Mole_frac
      real(kind=c_double), dimension(*), intent(out) :: Therm_factors
    end subroutine calc_therm_factor_loc
  end interface  

  !> This routine calculates Maxwell-Stefan diffusivity coeffcient Matrix
  !! for given mole_frac of all species
  interface
    subroutine calc_ms_diff_matrix_from_molefrac(nSpc, Temp, Press, Mole_frac, &
      & D_ij_out) bind(c, name='calc_ms_diff_matrix_from_molefrac_C')
      use, intrinsic :: iso_c_binding
      integer(kind=c_int), value, intent(in) :: nSpc
      real(kind=c_double), value, intent(in) :: Temp
      real(kind=c_double), value, intent(in) :: Press
      real(kind=c_double), dimension(*), intent(in) :: Mole_frac
      real(kind=c_double), dimension(*), intent(out) :: D_ij_out
    end subroutine calc_ms_diff_matrix_from_molefrac
  end interface

  !> This routine calculates Maxwell-Stefan diffusivity coeffcient Matrix
  !! for given mole_frac of all species
  interface
    subroutine calc_ms_diff_matrix_from_moledens(nSpc, Temp, Press, Mole_dens, &
      & D_ij_out) bind(c, name='calc_ms_diff_matrix_from_moledens_C')
      use, intrinsic :: iso_c_binding
      integer(kind=c_int), value, intent(in) :: nSpc
      real(kind=c_double), value, intent(in) :: Temp
      real(kind=c_double), value, intent(in) :: Press
      real(kind=c_double), dimension(*), intent(in) :: Mole_dens
      real(kind=c_double), dimension(*), intent(out) :: D_ij_out
    end subroutine calc_ms_diff_matrix_from_moledens
  end interface

  interface mus_calc_thermFactor
    module procedure mus_calc_thermFactor_single
  end interface mus_calc_thermFactor

  interface mus_calc_MS_DiffMatrix
    module procedure mus_calc_MS_DiffMatrix_single
  end interface

contains

  ! ************************************************************************** !
  !> This function loads property file using external c-function
  function mus_init_eNRTL( filename, nFields ) result(success)
    ! -------------------------------------------------------------------------!
    character(kind=c_char), dimension(*) :: filename
    !> number of fields in mixture
    integer, intent(out) :: nFields
    logical :: success
    ! -------------------------------------------------------------------------!
    integer(kind=c_int) :: nFields_c
    integer(kind=c_int) :: success_c
    ! -------------------------------------------------------------------------!

    success_c = init_eNRTL_loc( filename, nFields_c )
    if (success_c==0) then
      success = .true.
      nFields = nFields_c
    else
      success = .false.
      nFields = -1
    end if
     
  end function mus_init_eNRTL

  ! ************************************************************************** !
  !> This routine calculates thermodynamic factor for given mole_frac
  !! of all species for single element
  subroutine mus_calc_thermFactor_single( nFields, temp, press, mole_frac,     &
    &                                     therm_factors )
    ! -------------------------------------------------------------------------!
    !> number of fields in mixture
    integer, intent(in) :: nFields
    !> mixture temperature
    real(kind=rk), intent(in) :: temp
    !> mixture pressure
    real(kind=rk), intent(in) :: press
    !> mole fraction of all species of single element
    real(kind=rk), intent(in) :: mole_frac(nFields)
    !> thermodynamic factor matrix
    real(kind=rk), intent(out) :: therm_factors(nFields,nFields)
    ! -------------------------------------------------------------------------!
    integer(kind=c_int) :: nFields_c
    real(kind=c_double) :: temp_c, press_c
    real(kind=c_double) :: mole_frac_c(nFields)
    real(kind=c_double):: therm_factors_c(nFields*nFields)
    integer :: iField, iField_2
    ! -------------------------------------------------------------------------!

    nFields_c = int(nFields, kind=c_int)
    temp_c = real(temp, kind=c_double)
    press_c = real(press, kind=c_double)
    do iField = 1,nFields
      mole_frac_c(iField) = real(mole_frac(iField), kind=c_double)
    end do

    call calc_therm_factor_loc( nFields_c, temp_c, press_c, mole_frac_c,       &
      &                         therm_factors_c )
    do iField = 1, nFields
      do iField_2 = 1, nFields
        therm_factors(iField, iField_2) =                                      &
          &                    therm_factors_c((iField-1)*nFields+ iField_2)
      end do
    end do
!    therm_factors = 0.0_rk
!    do iField = 1, nFields
!      therm_factors(iField, iField) = 1.0_rk
!    end do

!    write(*,*) 'Therm_factors_c ', therm_factors_c  
!    write(*,*) 'Therm_factors '
!    do iField = 1, nFields
!    write(*,*) 'iField ', iField, ' factors ', therm_factors(iField,:)
!    end do
  end subroutine mus_calc_thermFactor_single  

  ! ************************************************************************** !
  !> This routine calculates Diffusivity coefficients matrix for given mole_frac
  !! of all species for single element
  subroutine mus_calc_MS_DiffMatrix_single( nFields, temp, press, mole_dens,   &
    &                                     D_ij_out )
    ! -------------------------------------------------------------------------!
    !> number of fields in mixture
    integer, intent(in) :: nFields
    !> mixture temperature
    real(kind=rk), intent(in) :: temp
    !> mixture pressure
    real(kind=rk), intent(in) :: press
    !> mole density of all species of single element
    real(kind=rk), intent(in) :: mole_dens(nFields)
    !> thermodynamic factor matrix
    real(kind=rk), intent(out) :: D_ij_out(nFields,nFields)
    ! -------------------------------------------------------------------------!
    integer(kind=c_int) :: nFields_c
    real(kind=c_double) :: temp_c, press_c
    real(kind=c_double) :: mole_dens_c(nFields)
    real(kind=c_double):: D_ij_out_c(nFields*nFields)
    integer :: iField, iField_2
    ! -------------------------------------------------------------------------!

    nFields_c = int(nFields, kind=c_int)
    temp_c = real(temp, kind=c_double)
    press_c = real(press, kind=c_double)
    do iField = 1,nFields
      mole_dens_c(iField) = real(mole_dens(iField), kind=c_double)
    end do

!    write(*,*) 'nFields ', nFields
!    write(*,*) 'temp ', temp
!    write(*,*) 'press ', press
!    write(*,*) 'mole_dens ', mole_dens
    
    call calc_ms_diff_matrix_from_moledens( nFields_c, temp_c, press_c,        &
      & mole_dens_c, D_ij_out_c )
    do iField = 1, nFields
      do iField_2 = 1, nFields
        D_ij_out(iField, iField_2) = D_ij_out_c((iField-1)*nFields+ iField_2)
      end do
    end do
!    write(*,*) 'D_ij_out_c ', D_ij_out_c  
!    write(*,*) 'D_ij_out '
!    do iField = 1, nFields
!      write(*,*) 'iField ', iField, ' factors ', D_ij_out(iField,:)
!    end do
  end subroutine mus_calc_MS_DiffMatrix_single

end module mus_eNRTL_module
