text:
!> Interface to thermodynamic models.
!! Using Helmholtz formulation in temperature and volume
!!
!! \author MH, 2015-02
module eosTV
  use thermopack_var, only: nc, nce, get_active_eos, base_eos_param, &
       thermo_model, get_active_thermo_model, apparent_to_real_mole_numbers, &
       real_to_apparent_diff
  use thermopack_constants
  !
  implicit none
  save
  !
  ! Include TREND interface
  include 'trend_interface.f95'
  !
  private
  public :: pressure, internal_energy_tv, free_energy_tv, entropy_tv
  public :: thermo_tv, Fres, Fideal, chemical_potential_tv
  public :: virial_coefficients, secondvirialcoeffmatrix
  public :: binaryThirdVirialCoeffMatrix
  public :: enthalpy_tv, Fres_ne, Fideal_ne
  !
  public :: enthalpy_tvp, entropy_tvp, thermo_tvp
contains

  !----------------------------------------------------------------------
  !> Calculate pressure given composition, temperature and density.
  !>
  !> \author MH, 2012-03-14
  !----------------------------------------------------------------------
  function pressure(t,v,n,dpdv,dpdt,d2pdv2,dpdn,recalculate) result(p)
    use single_phase, only: TV_CalcPressure
    use cbHelm, only: cbPvv, cbPi
    implicit none
    ! Transferred variables
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: v !< m3 - Volume
    real, dimension(1:nc), intent(in) :: n !< Mol numbers
    real, optional, intent(out) :: dpdt !< Pa/K - Pressure differential wrpt. temperature
    real, optional, intent(out) :: dpdv !< Pa/m3 - Pressure differential wrpt. specific volume
    real, optional, intent(out) :: d2pdv2 !< Pa/m6 - Second pressure differential wrpt. specific volume
    real, optional, dimension(nc), intent(out) :: dpdn !< Pa/mol - Second pressure differential wrpt. specific mole numbers
    logical, optional, intent(in) :: recalculate !< flag for Thermopack: if true, recalculate cbeos-structure
    real :: p !< Pa - Pressure
    ! Locals
    real :: d2Pdrho2, rho, dPdrho
    !integer :: VSHIFTID
    logical :: recalculate_loc
    class(base_eos_param), pointer :: act_eos_ptr
    type(thermo_model), pointer :: act_mod_ptr
    !--------------------------------------------------------------------
    if (present(recalculate)) then
      recalculate_loc = recalculate
    else
      recalculate_loc = .true.
    end if
    act_mod_ptr => get_active_thermo_model()
    select case (act_mod_ptr%EoSlib)
    case (THERMOPACK)
      ! Thermopack
      act_eos_ptr => get_active_eos()
      call TV_CalcPressure(nc,act_mod_ptr%comps,act_eos_ptr,T,v,n,p,&
           dpdv, dpdt, d2pdv2, dpdn, recalculate=recalculate_loc)
    case (TREND)
      ! TREND
      p = trend_pressure(n,t,v,dpdv,dpdt,dpdn)
      if (present(d2pdv2)) then
        rho = 1.0/v
        d2Pdrho2 = trend_d2Pdrho2(n,T,rho)
        if (present(dpdv)) then
          dPdrho = -dPdv*v**2
        else
          dPdrho = trend_dpdrho(n,T,rho)
        endif
        d2pdv2 = rho**4*d2Pdrho2 + 2.0*rho**3*dPdrho
      endif
    case default
      write(*,*) 'EoSlib error in eos::pressure: No such EoS libray:',act_mod_ptr%EoSlib
      call stoperror('')
    end select
  end function pressure

  !----------------------------------------------------------------------
  !> Calculate internal energy given composition, temperature and density.
  !>
  !> \author MH, 2012-03-14
  !----------------------------------------------------------------------
  subroutine internal_energy_tv(t,v,n,u,dudt,dudv,dudn,contribution)
    use eosdata
    use thermopack_constants, only: PROP_RESIDUAL, PROP_IDEAL
    implicit none
    ! Transferred variables
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: v !< m3 - Specific volume
    real, dimension(1:nc), intent(in) :: n !< Mol numbers
    real, intent(out) :: u !< J - Specific internal energy
    real, optional, intent(out) :: dudt !< J/K - Energy differential wrpt. temperature
    real, optional, intent(out) :: dudv !< J/m3 - Energy differential wrpt. volume
    real, optional, intent(out) :: dudn(nc) !< J/mol - Energy differential wrpt. mol numbers
    integer, optional, intent(in) :: contribution !< Contribution from ideal (PROP_IDEAL), residual (PROP_RESIDUAL) or both (PROP_OVERALL)
    ! Locals
    !--------------------------------------------------------------------
    real :: F_T, u_id, F_Tn(nc)
    real, pointer :: F_TT_p, F_TV_p, F_Tn_p(:)
    real, target :: F_TT_l, F_TV_l, F_Tn_l(nce)
    real, dimension(nce) :: ne
    logical :: res, ideal
    call apparent_to_real_mole_numbers(n,ne)

    res = .true.
    ideal = .true.
    if (present(contribution)) then
      if (contribution == PROP_RESIDUAL) then
        ideal = .false.
      else if (contribution == PROP_IDEAL) then
        res = .false.
      endif
    endif

    if (present(dudt)) then
      dudt = 0
      F_TT_p => F_TT_l
    else
      F_TT_p => NULL()
    endif
    if (present(dudv)) then
      dudv = 0
      F_TV_p => F_TV_l
    else
      F_TV_p => NULL()
    endif
    if (present(dudn)) then
      dudn = 0
      F_Tn_p => F_Tn_l
    else
      F_Tn_p => NULL()
    endif
    u = 0

    if (res) then
      call Fres_ne(T,V,ne,F_T=F_T,F_TT=F_TT_p,&
           F_TV=F_TV_p,F_Tn=F_Tn_p)
      u = (-Rgas)*T**2*F_T
      if (present(dudv)) dudv = (-Rgas)*T**2*F_TV_l
      if (present(dudt)) dudt = (-Rgas)*T*(T*F_TT_l+2.0*F_T)
      if (present(dudn)) then
        call real_to_apparent_diff(F_Tn_l,F_Tn)
        dudn = (-Rgas)*T**2*F_Tn
      endif
    endif
    if (ideal) then
      ! Add ideal gas contributions
      call Fideal_ne(T,V,ne,F_T=F_T,F_TT=F_TT_p,&
           F_TV=F_TV_p,F_Tn=F_Tn_p)
      u_id = (-Rgas)*T**2*F_T
      u = u + u_id
      if (present(dudv)) dudv = dudv + (-Rgas)*T**2*F_TV_l
      if (present(dudt)) dudt = dudt + (-Rgas)*T*(T*F_TT_l+2.0*F_T)
      if (present(dudn)) then
        call real_to_apparent_diff(F_Tn_l,F_Tn)
        dudn = dudn + (-Rgas)*T**2*F_Tn
      endif
    endif
  end subroutine internal_energy_tv

  !----------------------------------------------------------------------
  !> Calculate Helmholtz free energy given composition, temperature and density.
  !>
  !> \author GL, 2015-01-23
  !----------------------------------------------------------------------
  subroutine free_energy_tv(t,v,n,y,dydt,dydv,dydn,contribution)
    use thermopack_constants, only: Rgas, PROP_RESIDUAL, PROP_IDEAL
    implicit none
    ! Transferred variables
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: v !< m3 - Volume
    real, dimension(1:nc), intent(in) :: n !< Mol numbers
    real, intent(out) :: y !< J - Free energy
    real, optional, intent(out) :: dydt !< J/K - Differential wrt. temperature
    real, optional, intent(out) :: dydv !< J/m3 - Differential wrt. specific volume
    real, optional, intent(out) :: dydn(nc) !< J/mol - Helmholtz differential wrpt. mol numbers
    integer, optional, intent(in) :: contribution !< Contribution from ideal (PROP_IDEAL), residual (PROP_RESIDUAL) or both (PROP_OVERALL)
    ! Locals
    logical :: res, ideal
    real :: F, F_n(nc)
    real, pointer :: F_T_p, F_V_p, F_n_p(:)
    real, target :: F_T_l, F_V_l, F_n_l(nce)
    real, dimension(nce) :: ne
    !--------------------------------------------------------------------
    call apparent_to_real_mole_numbers(n,ne)

    res = .true.
    ideal = .true.
    if (present(contribution)) then
      if (contribution == PROP_RESIDUAL) then
        ideal = .false.
      else if (contribution == PROP_IDEAL) then
        res = .false.
      endif
    endif

    if (present(dYdT)) then
      dYdT = 0
      F_T_p => F_T_l
    else
      F_T_p => NULL()
    endif
    if (present(dYdV)) then
      dYdV = 0
      F_V_p => F_V_l
    else
      F_V_p => NULL()
    endif
    if (present(dYdn)) then
      dYdn = 0
      F_n_p => F_n_l
    else
      F_n_p => NULL()
    endif
    y = 0

    if (res) then
      call Fres_ne(T,V,ne,F=F,F_T=F_T_p,F_V=F_V_p,F_n=F_n_p)
      Y = Rgas*T*F
      if (present(dYdt)) dYdt = Rgas*(F + T*F_T_l)
      if (present(dYdv)) dYdv = Rgas*T*F_V_l
      if (present(dYdn)) then
        call real_to_apparent_diff(F_n_l,F_n)
        dYdn = Rgas*T*F_n
      endif
    endif
    if (ideal) then
      ! Add ideal gas contributions
      call Fideal_ne(T,V,ne,F=F,F_T=F_T_p,F_V=F_V_p,F_n=F_n_p)
      Y = Y + Rgas*T*F
      if (present(dYdt)) dYdt = dYdt + Rgas*(F + T*F_T_l)
      if (present(dYdv)) dYdv = dYdv + Rgas*T*F_V_l
      if (present(dYdn)) then
        call real_to_apparent_diff(F_n_l,F_n)
        dYdn = dYdn + Rgas*T*F_n
      endif
    endif
  end subroutine free_energy_tv

  !----------------------------------------------------------------------
  !> Calculate entropy given composition, temperature and density.
  !>
  !> \author MH, 2015-02
  !----------------------------------------------------------------------
  subroutine entropy_tv(t,v,n,s,dsdt,dsdv,dsdn,contribution)
    use thermopack_constants, only: Rgas, PROP_RESIDUAL, PROP_IDEAL
    implicit none
    ! Transferred variables
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: v !< m3 - Volume
    real, dimension(1:nc), intent(in) :: n !< Mol numbers
    real, intent(out) :: s !< J/K - Entropy
    real, optional, intent(out) :: dsdt !< J/K2 - Entropy differential wrpt. temperature
    real, optional, intent(out) :: dsdv !< J/K/m3 - Entropy differential wrpt. specific volume
    real, optional, intent(out) :: dsdn(nc) !< J/K/mol - Entropy differential wrpt. mol numbers
    integer, optional, intent(in) :: contribution !< Contribution from ideal (PROP_IDEAL), residual (PROP_RESIDUAL) or both (PROP_OVERALL)
    ! Locals
    real :: s_id, dsdt_id, dsdv_id, F, F_T, dsdne_id(nce), dsdne(nce)
    real :: ne(nce)
    logical :: res, ideal
    real, pointer :: F_TT_p, F_V_p, F_TV_p
    real, target :: F_TT, F_V, F_TV
    real, pointer :: F_Tne_p(:), F_ne_p(:)
    real, target :: F_Tne(nce), F_ne(nce)
    !--------------------------------------------------------------------
    call apparent_to_real_mole_numbers(n,ne)
    res = .true.
    ideal = .true.
    if (present(contribution)) then
      if (contribution == PROP_RESIDUAL) then
        ideal = .false.
      else if (contribution == PROP_IDEAL) then
        res = .false.
      endif
    endif

    if (present(dsdT)) then
      dsdT = 0
      F_TT_p => F_TT
    else
      F_TT_p => NULL()
    endif
    if (present(dsdV)) then
      dsdV = 0
      F_V_p => F_V
      F_TV_p => F_TV
    else
      F_V_p => NULL()
      F_TV_p => NULL()
    endif
    if (present(dsdn)) then
      dsdne = 0
      F_ne_p => F_ne
      F_Tne_p => F_Tne
    else
      F_ne_p => NULL()
      F_Tne_p => NULL()
    endif
    s = 0

    if (res) then
      ! Residual contribution
      call Fres_ne(T,v,ne,F=F,F_T=F_T,F_V=F_V_p,F_n=F_ne_p,F_TT=F_TT_p,&
           F_TV=F_TV_p,F_Tn=F_Tne_p)
      s = - Rgas*(F + T*F_T)
      if (present(dsdt)) then
        dsdt = - Rgas*(2*F_T + T*F_TT)
      endif
      if (present(dsdv)) then
        dsdv = - Rgas*(F_V + T*F_TV)
      endif
      if(present(dsdn)) then
        dsdne = -Rgas*(F_ne + T*F_Tne)
      end if
    endif

    if (ideal) then
      call Fideal_ne(T,v,ne,F,F_T,F_V=F_V_p,F_n=F_ne_p,F_TT=F_TT_p,&
           F_TV=F_TV_p,F_Tn=F_Tne_p)
      s_id = - Rgas*(F + T*F_T)
      s = s + s_id
      if (present(dsdt)) then
        dsdt_id = - Rgas*(2*F_T + T*F_TT)
        dsdt = dsdt + dsdt_id
      endif
      if (present(dsdv)) then
        dsdv_id = - Rgas*(F_V + T*F_TV)
        dsdv = dsdv + dsdv_id
      endif
      if(present(dsdn)) then
        dsdne_id = -Rgas*(F_ne + T*F_Tne)
        dsdne = dsdne + dsdne_id
      end if
    endif
    if (present(dSdn)) then
      call real_to_apparent_diff(dSdne,dSdn)
    endif
  end subroutine entropy_tv

  !----------------------------------------------------------------------
  !> Calculate enthalpy given composition, temperature and density.
  !>
  !> \author MH, 2019-06
  !----------------------------------------------------------------------
  subroutine enthalpy_tv(t,v,n,h,dhdt,dhdv,dhdn,contribution)
    use thermopack_constants, only: Rgas
    implicit none
    ! Transferred variables
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: v !< m3 - Volume
    real, dimension(1:nc), intent(in) :: n !< Mol numbers
    real, intent(out) :: h !< J - Enthalpy
    real, optional, intent(out) :: dhdt !< J/K - Enthalpy differential wrpt. temperature
    real, optional, intent(out) :: dhdv !< J/m3 - Enthalpy differential wrpt. volume
    real, optional, intent(out) :: dhdn(nc) !< J/m3 - Enthalpy differential wrpt. mol numbers
    integer, optional, intent(in) :: contribution !< Contribution from ideal (PROP_IDEAL), residual (PROP_RESIDUAL) or both (PROP_OVERALL)
    ! Locals
    real :: h_id, dhdt_id, dhdv_id, F, F_T, F_V, h_T, h_id_T, dhdne_id(nc), dhdne(nce)
    real :: ne(nce)
    logical :: res, ideal
    real, pointer :: F_TT_p, F_TV_p, F_VV_p
    real, target :: F_TT, F_TV, F_VV
    real, pointer :: F_Tne_p(:), F_Vne_p(:)
    real, target :: F_Tne(nce), F_Vne(nce)
    !--------------------------------------------------------------------
    call apparent_to_real_mole_numbers(n,ne)
    res = .true.
    ideal = .true.
    if (present(contribution)) then
      if (contribution == PROP_RESIDUAL) then
        ideal = .false.
      else if (contribution == PROP_IDEAL) then
        res = .false.
      endif
    endif

    if (present(dhdT)) then
      dhdT = 0
      F_TT_p => F_TT
    else
      F_TT_p => NULL()
    endif
    if (present(dhdV)) then
      dhdV = 0
      F_VV_p => F_VV
      F_TV_p => F_TV
    else
      F_VV_p => NULL()
      F_TV_p => NULL()
    endif
    if (present(dhdn)) then
      dhdne = 0
      F_Vne_p => F_Vne
      F_Tne_p => F_Tne
    else
      F_Vne_p => NULL()
      F_Tne_p => NULL()
    endif
    h = 0

    if (res) then
      ! Residual contribution
      call Fres_ne(T,v,ne,F=F,F_T=F_T,F_V=F_V,F_TT=F_TT_p,&
           F_TV=F_TV_p,F_VV=F_VV_p,F_Tn=F_Tne_p,F_Vn=F_Vne_p)
      h_T = -Rgas*(T*F_T + v*F_V)
      h = h_T*T
      if (present(dhdt)) then
        dhdt = h_T - Rgas*T*(F_T + T*F_TT + v*F_TV)
      endif
      if (present(dhdv)) then
        dhdv = - Rgas*T*(T*F_TV + F_V + v*F_VV)
      endif
      if (present(dhdn)) then
        dhdne = - Rgas*T*(T*F_Tne + v*F_Vne)
      endif
    endif
    if (ideal) then
      call Fideal_ne(T,v,ne,F,F_T,F_V=F_V,F_TT=F_TT_p,&
           F_TV=F_TV_p,F_VV=F_VV_p,F_Tn=F_Tne_p,F_Vn=F_Vne_p)
      h_id_T = -Rgas*(T*F_T + v*F_V)
      h_id = h_id_T*T
      h = h + h_id
      if (present(dhdt)) then
        dhdt_id = h_id_T - Rgas*T*(F_T + T*F_TT + v*F_TV)
        dhdt = dhdt + dhdt_id
      endif
      if (present(dhdv)) then
        dhdv_id = - Rgas*T*(T*F_TV + F_V + v*F_VV)
        dhdv = dhdv + dhdv_id
      endif
      if (present(dhdn)) then
        dhdne_id = - Rgas*T*(T*F_Tne + v*F_Vne)
        dhdne = dhdne + dhdne_id
      endif
    endif
    if (present(dhdn)) then
      call real_to_apparent_diff(dhdne,dhdn)
    endif
  end subroutine enthalpy_tv

  !----------------------------------------------------------------------
  !> Calculate fugacity and differentials given composition,
  !> temperature and specific volume
  !>
  !> \author MH, 2015-10
  !----------------------------------------------------------------------
  subroutine thermo_tv(t,v,n,lnphi,lnphit,lnphiv,lnphin)
    use single_phase, only: TV_CalcFugacity
    !$ use omp_lib
    implicit none
    ! Transferred variables
    real, intent(in)                                  :: t !< K - Temperature
    real, intent(in)                                  :: v !< m3 - Volume
    real, dimension(1:nc), intent(in)                 :: n !< Mole numbers [mol]
    real, dimension(1:nc), intent(out)                :: lnphi !< Logarithm of fugasity
    real, optional, dimension(1:nc), intent(out)      :: lnphit !< 1/K - Logarithm of fugasity differential wrpt. temperature
    real, optional, dimension(1:nc), intent(out)      :: lnphiv !< mol/m3 - Logarithm of fugasity differential wrpt. volume
    real, optional, dimension(1:nc,1:nc), intent(out) :: lnphin !< Logarithm of fugasity differential wrpt. mole numbers
    ! Locals
    class(base_eos_param), pointer :: act_eos_ptr
    type(thermo_model), pointer :: act_mod_ptr
    !
    act_mod_ptr => get_active_thermo_model()
    select case (act_mod_ptr%EoSlib)
    case (THERMOPACK)
      ! Thermopack
      act_eos_ptr => get_active_eos()
      call TV_CalcFugacity(nc,act_mod_ptr%comps,act_eos_ptr,T,v,n,lnphi,&
           lnphiT,lnphiV,lnphin)
    case (TREND)
      ! TREND
      call trend_thermoTV(T,v,n,lnphi,lnphiT,lnphiV,lnphin)
    case default
      write(*,*) 'EoSlib error in eosTV::thermo: No such EoS libray:',act_mod_ptr%EoSlib
      call stoperror('')
    end select
  end subroutine thermo_tv

  !----------------------------------------------------------------------
  !> Calculate residual reduced Helmholtz energy
  !>
  !> \author MH, 2012-01-27
  !----------------------------------------------------------------------
  subroutine Fres(T,V,n,F,F_T,F_V,F_n,F_TT,F_TV,F_VV,F_Tn,F_Vn,F_nn,F_VVV,recalculate)
    use single_phase, only: TV_CalcFres
    use thermopack_var, only: nce, apparent_to_real_mole_numbers, real_to_apparent_differentials
    implicit none
    ! Transferred variables
    real, intent(in)                                  :: t !< K - Temperature
    real, intent(in)                                  :: v !< m3 - Volume
    real, dimension(1:nc), intent(in)                 :: n !< mol numbers
    real, optional, intent(out)                       :: F,F_T,F_V,F_TT,F_TV,F_VV,F_VVV !<
    real, optional, dimension(1:nc), intent(out)      :: F_n,F_Tn,F_Vn !<
    real, optional, dimension(1:nc,1:nc), intent(out) :: F_nn !<
    logical, optional, intent(in)                     :: recalculate
    ! Locals
    real, dimension(nce) :: ne
    real, target, dimension(nce) :: F_ne, F_Tne, F_Vne
    real, target, dimension(nce,nce) :: F_nene
    real, pointer :: F_ne_p(:), F_Tne_p(:), F_Vne_p(:)
    real, pointer :: F_nene_p(:,:)
    !
    !--------------------------------------------------------------------
    call apparent_to_real_mole_numbers(n,ne)
    if (present(F_n)) then
      F_ne_p => F_ne
    else
      F_ne_p => NULL()
    endif
    if (present(F_Vn)) then
      F_Vne_p => F_Vne
    else
      F_Vne_p => NULL()
    endif
    if (present(F_Tn)) then
      F_Tne_p => F_Tne
    else
      F_Tne_p => NULL()
    endif
    if (present(F_nn)) then
      F_nene_p => F_nene
    else
      F_nene_p => NULL()
    endif
    call Fres_ne(T,V,ne,F,F_T=F_T,F_V=F_V,F_n=F_ne_p,F_TT=F_TT,F_TV=F_TV,F_VV=F_VV,&
         F_Tn=F_Tne_p,F_Vn=F_Vne_p,F_nn=F_nene_p,F_VVV=F_VVV,recalculate=recalculate)
    call real_to_apparent_differentials(F_ne,F_Tne,F_Vne,F_nene,&
           F_n,F_Tn,F_Vn,F_nn)
  end subroutine Fres

  !----------------------------------------------------------------------
  !> Calculate residual reduced Helmholtz energy for real (not apparent) mol numbers
  !>
  !> \author MH, 2022-02
  !----------------------------------------------------------------------
  subroutine Fres_ne(T,V,ne,F,F_T,F_V,F_n,F_TT,F_TV,F_VV,F_Tn,F_Vn,F_nn,F_VVV,recalculate)
    use single_phase, only: TV_CalcFres
    use thermopack_var, only: nce, apparent_to_real_mole_numbers, real_to_apparent_differentials
    implicit none
    ! Transferred variables
    real, intent(in)                                    :: t !< K - Temperature
    real, intent(in)                                    :: v !< m3 - Volume
    real, dimension(1:nce), intent(in)                  :: ne !< mol numbers
    real, optional, intent(out)                         :: F,F_T,F_V,F_TT,F_TV,F_VV,F_VVV !<
    real, optional, dimension(1:nce), intent(out)       :: F_n,F_Tn,F_Vn !<
    real, optional, dimension(1:nce,1:nce), intent(out) :: F_nn !<
    logical, optional, intent(in)                       :: recalculate
    ! Locals
    class(base_eos_param), pointer :: act_eos_ptr
    type(thermo_model), pointer :: act_mod_ptr
    !
    !--------------------------------------------------------------------
    act_mod_ptr => get_active_thermo_model()
    select case (act_mod_ptr%EoSlib)
    case (THERMOPACK)
      ! Thermopack
      act_eos_ptr => get_active_eos()
      call TV_CalcFres(nce,act_mod_ptr%comps,act_eos_ptr,&
           T,V,ne,F=F,F_T=F_T,F_V=F_V,F_n=F_n,&
           F_TT=F_TT,F_TV=F_TV,F_VV=F_VV,F_Tn=F_Tn,F_Vn=F_Vn,F_nn=F_nn,&
           F_VVV=F_VVV,recalculate=recalculate)
    case (TREND)
      ! TREND
      call trend_calcFres(T,v,ne,F,F_T,F_V,F_n,F_TT,F_TV,F_VV,F_Tn,F_Vn,F_nn)
    case default
      write(*,*) 'EoSlib error in eosTV::Fres: No such EoS libray:',act_mod_ptr%EoSlib
      call stoperror('')
    end select
  end subroutine Fres_ne

  !----------------------------------------------------------------------
  !> Calculate ideal reduced Helmholtz energy
  !>
  !> \author MH, 2012-01-27
  !----------------------------------------------------------------------
  subroutine Fideal(T,V,n,F,F_T,F_V,F_n,F_TT,F_TV,F_VV,F_Tn,F_Vn,F_nn)
    use single_phase, only: TV_CalcFid
    use thermopack_var, only: nce, apparent_to_real_mole_numbers, real_to_apparent_differentials
    implicit none
    ! Transferred variables
    real, intent(in)                                  :: t !< K - Temperature
    real, intent(in)                                  :: v !< m3/mol - Volume
    real, dimension(1:nc), intent(in)                 :: n !< Compozition
    real, optional, intent(out)                       :: F,F_T,F_V,F_TT,F_TV,F_VV !<
    real, optional, dimension(1:nc), intent(out)      :: F_n,F_Tn,F_Vn !<
    real, optional, dimension(1:nc,1:nc), intent(out) :: F_nn !<
    ! Locals
    real, dimension(nce) :: ne
    real, target, dimension(nce) :: F_ne, F_Tne, F_Vne
    real, target, dimension(nce,nce) :: F_nene
    real, pointer :: F_ne_p(:), F_Tne_p(:), F_Vne_p(:)
    real, pointer :: F_nene_p(:,:)
    !
    !--------------------------------------------------------------------
    call apparent_to_real_mole_numbers(n,ne)
    if (present(F_n)) then
      F_ne_p => F_ne
    else
      F_ne_p => NULL()
    endif
    if (present(F_Vn)) then
      F_Vne_p => F_Vne
    else
      F_Vne_p => NULL()
    endif
    if (present(F_Tn)) then
      F_Tne_p => F_Tne
    else
      F_Tne_p => NULL()
    endif
    if (present(F_nn)) then
      F_nene_p => F_nene
    else
      F_nene_p => NULL()
    endif
    call Fideal_ne(T,V,ne,F=F,F_T=F_T,F_V=F_V,F_n=F_ne,&
         F_TT=F_TT,F_TV=F_TV,F_VV=F_VV,F_Tn=F_Tne_p,F_Vn=F_Vne_p,F_nn=F_nene_p)
    call real_to_apparent_differentials(F_ne,F_Tne,F_Vne,F_nene,&
         F_n,F_Tn,F_Vn,F_nn)
  end subroutine Fideal

  !----------------------------------------------------------------------
  !> Calculate ideal reduced Helmholtz energy for real (not apparent) components
  !>
  !> \author MH, 2022-02
  !----------------------------------------------------------------------
  subroutine Fideal_ne(T,V,ne,F,F_T,F_V,F_n,F_TT,F_TV,F_VV,F_Tn,F_Vn,F_nn)
    use single_phase, only: TV_CalcFid
    use thermopack_var, only: nce, apparent_to_real_mole_numbers, real_to_apparent_differentials
    implicit none
    ! Transferred variables
    real, intent(in)                                    :: t !< K - Temperature
    real, intent(in)                                    :: v !< m3/mol - Volume
    real, dimension(1:nce), intent(in)                  :: ne !< Compozition
    real, optional, intent(out)                         :: F,F_T,F_V,F_TT,F_TV,F_VV !<
    real, optional, dimension(1:nce), intent(out)       :: F_n,F_Tn,F_Vn !<
    real, optional, dimension(1:nce,1:nce), intent(out) :: F_nn !<
    ! Locals
    class(base_eos_param), pointer :: act_eos_ptr
    type(thermo_model), pointer :: act_mod_ptr
    !
    !--------------------------------------------------------------------
    act_mod_ptr => get_active_thermo_model()
    select case (act_mod_ptr%EoSlib)
    case (THERMOPACK)
      ! Thermopack
      act_eos_ptr => get_active_eos()
      call TV_CalcFid(nce,act_mod_ptr%comps,act_eos_ptr,T,V,ne,F=F,F_T=F_T,F_V=F_V,F_n=F_n,&
           F_TT=F_TT,F_TV=F_TV,F_VV=F_VV,F_Tn=F_Tn,F_Vn=F_Vn,F_nn=F_nn)
    case (TREND)
      ! TREND
      call trend_CalcFid(T,V,ne,F,F_T,F_V,F_TT,F_TV,F_VV,F_n,F_Tn,F_Vn,F_nn)
    case default
      write(*,*) 'EoSlib error in eosTV::Fideal: No such EoS libray:',act_mod_ptr%EoSlib
      call stoperror('')
    end select
  end subroutine Fideal_ne

  !----------------------------------------------------------------------
  !> Calculate (composition-dependent) virial coefficients B and C,
  !> defined as P/RT = rho + B*rho**2 + C*rho**3 + O(rho**4) as rho->0.
  !----------------------------------------------------------------------
  subroutine virial_coefficients(T,n,B,C)
    use single_phase, only: TV_CalcFres
    implicit none
    ! Transferred variables
    real, intent(in)                      :: t !< Temperature [K]
    real, dimension(1:nc), intent(in)     :: n !< Composition [-]
    real, intent(out)                     :: B !< Second virial coefficients [m3/mol]
    real, intent(out)                     :: C !< Third virial coefficient [m6/mol2]
    ! Locals
    real :: V, F_V, F_VV, z(nc)
    !
    !-------------------------------------------------------------------
    V = 1.0e7 ! might need tuning
    z = n/sum(n)
    call Fres(T,V,Z,F_V=F_V,F_VV=F_VV)
    B = -V**2*F_V
    C = V**4*(F_VV + 2.0*F_V/V)
  end subroutine virial_coefficients

  !----------------------------------------------------------------------
  !> Calculate composition-independent virial coefficients B,
  !! defined as P = RT*rho + B*rho**2 + C*rho**3 + O(rho**4) as rho->0.
  !! Including cross coefficients.
  !----------------------------------------------------------------------
  subroutine secondVirialCoeffMatrix(T,Bmat)
    implicit none
    ! Transferred variables
    real, intent(in)                      :: t !< Temperature [K]
    real, intent(out)                     :: Bmat(nc,nc) !< Second virial coefficients [m3/mol]
    ! Locals
    real :: Bii, Bij, C, z(nc), x
    integer :: i, j
    !
    !-------------------------------------------------------------------

    ! Pure virial coefficients
    do i=1,nc
      z = 0
      z(i) = 1.0
      call virial_coefficients(T,z,Bii,C)
      Bmat(i,i) = Bii
    end do

    ! Cross virial coefficients
    x = 0.5 ! The answer should be independent of the value of x
    do i=1,nc
      do j=i+1,nc
        z = 0
        z(i) = x
        z(j) = 1-x
        call virial_coefficients(T,z,Bij,C)
        Bmat(i,j) = (Bij - Bmat(i,i)*x**2 - Bmat(j,j)*(1-x)**2)/(2*x*(1-x))
        Bmat(j,i) = Bmat(i,j)
      end do
    end do
  end subroutine secondVirialCoeffMatrix

  !----------------------------------------------------------------------
  !> Calculate composition-independent virial coefficients C,
  !! defined as P = RT*rho + B*rho**2 + C*rho**3 + O(rho**4) as rho->0.
  !! Including cross coefficients
  !! Currently the code only support binary mixtures
  !----------------------------------------------------------------------
  subroutine binaryThirdVirialCoeffMatrix(T,Cmat)
    implicit none
    ! Transferred variables
    real, intent(in)                      :: t !< Temperature [K]
    real, intent(out)                     :: Cmat(nc,nc) !< Third virial coefficients [m6/mol2]
    ! Locals
    real :: Ciii, Cij1, Cij2, B, z(nc), x
    real :: A(2,2), D(2)
    integer :: i, j
    integer :: ifail, ipiv(2)
    logical :: lufail
    !
    !-------------------------------------------------------------------
    if (nc /= 2) then
      call stoperror("binaryThirdVirialCoeffMatrix only intended binary systems")
    endif
    ! Pure virial coefficients
    do i=1,nc
      z = 0
      z(i) = 1.0
      call virial_coefficients(T,z,B,Ciii)
      Cmat(i,i) = Ciii
    end do

    ! Cross virial coefficients for binaries
    x = 0.75
    do i=1,nc
      do j=i+1,nc
        z = 0
        z(i) = x
        z(j) = 1-x
        call virial_coefficients(T,z,B,Cij1)
        A(1,1) = 3.0*z(i)**2*z(j)
        A(1,2) = 3.0*z(j)**2*z(i)
        D(1) = -Cij1 + Cmat(i,i)*z(i)**3 + Cmat(j,j)*z(j)**3
        z(j) = x
        z(i) = 1-x
        call virial_coefficients(T,z,B,Cij2)
        A(2,1) = 3.0*z(i)**2*z(j)
        A(2,2) = 3.0*(1-x)*z(j)**2*z(i)
        D(2) = -Cij2 + Cmat(i,i)*z(i)**3 + Cmat(j,j)*z(j)**3

        ! Solve linear system
        lufail = .false.
        call dgetrf(2,2,A,2,ipiv,ifail)
        if (ifail /= 0) lufail = .true.
        if (.not. lufail) then
          call dgetrs('n',2,1,A,2,ipiv,D,2,ifail)
          if (ifail /= 0) lufail = .true.
        endif
        Cmat(i,j) = D(1)
        Cmat(j,i) = D(2)
      end do
    end do

  end subroutine binaryThirdVirialCoeffMatrix

  !> Calculate chemical potential and derivatives
  !>
  !> \author MAG, 2018-10-31
  !----------------------------------------------------------------------
  subroutine chemical_potential_tv(t, v, n, mu, dmudt, dmudv, dmudn, contribution)
    use thermopack_constants, only: rgas
    use thermopack_var, only: real_to_apparent_differentials
    implicit none
    real,                             intent(in)  :: t !< K - Temperature
    real,                             intent(in)  :: v !< m3 - Molar volume
    real, dimension(nc),              intent(in)  :: n !< mol - Mol numbers
    real, dimension(nc),              intent(out) :: mu !< J/mol
    real, dimension(nc),    optional, intent(out) :: dmudv !< J/m^3
    real, dimension(nc),    optional, intent(out) :: dmudt !< J/mol K
    real, dimension(nc,nc), optional, intent(out) :: dmudn !< J/mol^2
    integer, optional, intent(in) :: contribution !< Contribution from ideal (PROP_IDEAL), residual (PROP_RESIDUAL) or both (PROP_OVERALL)
    ! Locals
    logical :: res, ideal
    real :: Fr_ne(nce), Fi_ne(nce), mue(nce), dmudt_e(nce), dmudv_e(nce)
    real :: ne(nce), dmudne_e(nce,nce)
    real, pointer :: Fr_Tne_p(:), Fr_vne_p(:), Fr_nene_p(:,:)
    real, target :: Fr_Tne(nce), Fr_vne(nce), Fr_nene(nc,nce)
    real, pointer :: Fi_Tne_p(:), Fi_vne_p(:), Fi_nene_p(:,:)
    real, target :: Fi_Tne(nce), Fi_vne(nce), Fi_nene(nce,nce)
    !
    call apparent_to_real_mole_numbers(n,ne)
    res = .true.
    ideal = .true.
    if (present(contribution)) then
      if (contribution == PROP_RESIDUAL) then
        ideal = .false.
      else if (contribution == PROP_IDEAL) then
        res = .false.
      endif
    endif

    if (present(dmudt)) then
      dmudt_e = 0
      Fr_Tne_p => Fr_Tne
      Fi_Tne_p => Fi_Tne
    else
      Fr_Tne_p => NULL()
      Fi_Tne_p => NULL()
    endif
    if (present(dmudV)) then
      dmudV_e = 0
      Fr_Vne_p => Fr_Vne
      Fi_Vne_p => Fi_Vne
    else
      Fr_Vne_p => NULL()
      Fi_Vne_p => NULL()
    endif
    if (present(dmudn)) then
      dmudne_e = 0
      Fr_nene_p => Fr_nene
      Fi_nene_p => Fi_nene
    else
      Fr_nene_p => NULL()
      Fi_nene_p => NULL()
    endif
    mue = 0

    if (res) then
      call Fres_ne(t, v, ne, f_n=Fr_ne, F_Tn=Fr_Tne_p, F_Vn=Fr_Vne_p, F_nn=Fr_nene_p)
      mue = rgas*t*Fr_ne
      !
      if (present(dmudt)) then
        dmudt_e = rgas*(Fr_ne + t*Fr_Tne)
      end if
      if (present(dmudv)) then
        dmudv_e = rgas*t*Fr_Vne
      end if
      if (present(dmudn)) then
        dmudne_e = rgas*t*Fr_nene
      end if
    endif
    if (ideal) then
      call Fideal_ne(t, v, ne, f_n=Fi_ne, F_Tn=Fi_Tne_p, F_Vn=Fi_Vne_p, F_nn=Fi_nene_p)
      mue = mue + rgas*t*Fi_ne
      !
      if (present(dmudt)) then
        dmudt_e = dmudt_e + rgas*(Fi_ne + t*Fi_Tne)
      end if
      if (present(dmudv)) then
        dmudv_e = dmudv_e + rgas*t*Fi_Vne
      end if
      if (present(dmudn)) then
        dmudne_e = dmudne_e + rgas*t*Fi_nene
      end if
    endif
    !
    call real_to_apparent_differentials(mue,dmudt_e,dmudv_e,dmudne_e,&
         mu,dmudt,dmudv,dmudn)
  end subroutine chemical_potential_tv

  !-----------------------------------------------------------------------------
  !> Calculate the logarithmic fugacity coefficient and its differentials.
  !! Evaluate using (T,v,n) but output differentials for  (T, P, n)
  !!
  !! \author Morten Hammer, 2022-02
  !-----------------------------------------------------------------------------
  subroutine thermo_tvp(T,v,n,lnfug,dlnfugdT,dlnfugdP,dlnfugdn)
    use thermopack_constants, only: Rgas
    use thermopack_var, only: nce, apparent_to_real_mole_numbers, &
         TP_lnfug_apparent, base_eos_param
    ! Input.
    real, intent(in) :: T                               !< Temperature [K]
    real, intent(in) :: v                               !< Volume [m3]
    real, intent(in) :: n(nc)                           !< Mole numbers [mols]
    ! Output
    real, intent(out) :: lnfug(nc)
    real, optional, intent(out) :: dlnfugdt(nc), dlnfugdp(nc), dlnfugdn(nc,nc)
    ! Locals.
    real :: P     !< Pressure [Pa].
    real :: sumne  !< Total mole number in mixture [mol]
    real :: zFac
    real :: dPdV, dPdT, dPdn(nce)
    real :: dVdn(nce), F_n(nce), F_V
    real :: ne(nce), lnfug_real(nce)
    real :: dlnfugdt_real(nce), dlnfugdp_real(nce), dlnfugdn_real(nce,nce)
    integer :: i,j
    logical :: differentials
    real, target :: F_Tn(nce),F_TV,F_VV,F_Vn(nce),F_nn(nce,nce)
    real, pointer :: F_Tn_p(:),F_TV_p,F_VV_p,F_Vn_p(:),F_nn_p(:,:)

    call apparent_to_real_mole_numbers(n,ne)
    sumne = sum(ne)

    F_TV_p => NULL()
    F_VV_p => NULL()
    F_nn_p => NULL()
    F_Tn_p => NULL()
    F_Vn_p => NULL()
    differentials = .false.
    if (present(dlnfugdt) .or. present(dlnfugdp) .or. present(dlnfugdn)) then
      differentials = .true.
      F_Vn_p => F_Vn
      F_VV_p => F_VV
      if (present(dlnfugdt)) then
        F_Tn_p => F_Tn
        F_TV_p => F_TV
      endif
      if (present(dlnfugdn)) then
        F_nn_p => F_nn
      endif
    endif

    call Fres_ne(T=T,V=V,ne=ne,F_n=F_n,F_V=F_V,&
         F_VV=F_VV,F_Vn=F_Vn,F_TV=F_TV,F_Tn=F_Tn,F_nn=F_nn)
    P = -Rgas*T*F_V + sumne*Rgas*T/V
    zFac = V*P/(sumne*Rgas*T)
    if (differentials) then
      dPdV = -Rgas*T*(F_VV + sumne/V**2)
      dPdn = Rgas*T*(-F_Vn + 1/V)
      dVdn = -dPdn/dPdV
    end if

    lnfug_real = F_n - log(zFac)

    if (present(dlnfugdt)) then
      dPdT = P/T-Rgas*T*F_TV
      dlnfugdt_real = F_Tn + (1 - dVdn*dPdT/Rgas)/T
    endif

    if (present(dlnfugdp)) then
      dlnfugdp_real = dVdn/(Rgas*T)-1/P
    endif

    if (present(dlnfugdn)) then
      do i=1,nce
        do j=1,nce
          dlnfugdn_real(i,j) = F_nn(i,j) + 1/sumne - dVdn(j)*dPdn(i)/(Rgas*T)
        end do
      end do
    endif

    call TP_lnfug_apparent(nc,ne,n,P,lnfug_real,lnfug,dlnfugdt_real,&
       dlnfugdp_real,dlnfugdn_real,dlnfugdT,dlnfugdP,dlnfugdn)

  end subroutine thermo_tvp

  !----------------------------------------------------------------------
  !> Calculate enthalpy given composition, temperature and density.
  !> Differentials at constanr pressure
  !> \author MH, 2019-06
  !----------------------------------------------------------------------
  subroutine enthalpy_tvp(t,v,n,h,dhdt,dhdp,dhdn,contribution)
    use thermopack_constants, only: Rgas, PROP_RESIDUAL, PROP_IDEAL
    use ideal, only: idealEnthalpySingle
    use thermopack_var, only: nce, apparent_to_real_mole_numbers, &
         real_to_apparent_diff
    implicit none
    ! Transferred variables
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: v !< m3 - Volume
    real, dimension(1:nc), intent(in) :: n !< Mol numbers
    real, intent(out) :: h !< J - Enthalpy
    real, optional, intent(out) :: dhdt !< J/K - Enthalpy differential wrpt. temperature  (const pressure)
    real, optional, intent(out) :: dhdp !< J/Pa - Enthalpy differential wrpt. pressure
    real, optional, intent(out) :: dhdn(nc) !< J/m3 - Enthalpy differential wrpt. mol numbers
    integer, optional, intent(in) :: contribution !< Contribution from ideal (PROP_IDEAL), residual (PROP_RESIDUAL) or both (PROP_OVERALL)
    ! Locals
    real :: h_id, dhdt_id, F, F_T, F_V, h_T, dhdn_id(nc)
    real :: sumne, ne(nce), P, dPdV, dPdT, dVdT, dPdn(nce), dVdn(nce), dhdne(nce)
    integer :: j
    logical :: res, ideal_gas, differentials
    real, pointer :: F_TT_p, F_TV_p, F_VV_p
    real, target :: F_TT, F_TV, F_VV
    real, pointer :: F_Tn_p(:), F_Vn_p(:), F_n_p(:)
    real, target :: F_Tn(nce), F_Vn(nce), F_n(nce)
    !--------------------------------------------------------------------
    res = .true.
    ideal_gas = .true.
    if (present(contribution)) then
      if (contribution == PROP_RESIDUAL) then
        ideal_gas = .false.
      else if (contribution == PROP_IDEAL) then
        res = .false.
      endif
    endif

    F_TV_p => NULL()
    F_VV_p => NULL()
    F_TT_p => NULL()
    F_n_p => NULL()
    F_Tn_p => NULL()
    F_Vn_p => NULL()
    differentials = .false.
    h = 0
    if (present(dhdT) .or. present(dhdp) .or. present(dhdn)) then
      differentials = .true.
      F_TV_p => F_TV
      F_VV_p => F_VV
      if (present(dhdT)) then
        F_TT_p => F_TT
        dHdt = 0
      endif
      if (present(dhdP)) then
        dHdp = 0
      endif
      if (present(dhdn)) then
        F_n_p => F_n
        F_Tn_p => F_Tn
        F_Vn_p => F_Vn
        dHdn = 0
      endif
    endif

    ! Get real mol numbers
    call apparent_to_real_mole_numbers(n,ne)
    sumne = sum(ne)

    if (res) then
      ! Residual contribution
      call Fres_ne(T,v,ne,F=F,F_T=F_T,F_V=F_V,F_TT=F_TT_p,&
           F_TV=F_TV_p,F_VV=F_VV_p,F_Tn=F_Tn_p,F_Vn=F_Vn_p)
      h_T = -Rgas*(T*F_T + v*F_V)
      h = h_T*T
      if (differentials) then
        P = -Rgas*T*F_V + sumne*Rgas*T/V
        dPdV = -Rgas*T*(F_VV + sumne/V**2)
        dPdT = P/T-Rgas*T*F_TV
        dVdT = -dPdT/dPdV
        if (present(dHdt)) then
          dHdt = T*(dVdT*dPdT - Rgas*(2*F_T + T*F_TT + sumne/T))
        endif
        if (present(dHdp)) then
          dHdp = V-T*dVdt
        endif
        if (present(dHdn)) then
          dPdn = Rgas*T*(-F_Vn + 1/V)
          dVdn = -dPdn/dPdV
          dHdne = T*(dVdn*dPdt - Rgas*(F_Tn*T + 1))
          call real_to_apparent_diff(dhdne,dhdn)
        endif
      endif
    endif

    if (ideal_gas) then
      dhdt_id = 0
      do j=1,nc
        call idealEnthalpySingle(t,j,h_id,F_TT_p)
        h = h + n(j)*h_id
        if (present(dhdt)) then
          dhdt_id = dhdt_id + n(j)*F_TT
        endif
        ! No contribution to dhdp
        if (present(dhdn)) then
          dhdn_id(j) = h_id
        endif
      enddo
      if (present(dhdt)) then
        dhdt = dhdt + dhdt_id
      endif
      if (present(dhdn)) then
        dhdn = dhdn + dhdn_id
      endif
    endif

  end subroutine enthalpy_tvp

  !----------------------------------------------------------------------
  !> Calculate entropy given composition, temperature and density.
  !> Differentials at constant pressure
  !> \author MH, 2022-02
  !----------------------------------------------------------------------
  subroutine entropy_tvp(t,v,n,s,dsdt,dsdp,dsdn,contribution)
    use thermopack_constants, only: Rgas, PROP_RESIDUAL, PROP_IDEAL
    use ideal, only: idealEntropy_ne
    use thermopack_var, only: nce, apparent_to_real_mole_numbers, &
         real_to_apparent_diff
    implicit none
    ! Transferred variables
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: v !< m3 - Volume
    real, dimension(1:nc), intent(in) :: n !< Mol numbers
    real, intent(out) :: s !< J/K - Entropy
    real, optional, intent(out) :: dsdt !< J/K2 - Entropy differential wrpt. temperature (const pressure)
    real, optional, intent(out) :: dsdp !< J/K/Pa - Entropy differential wrpt. specific pressure
    real, optional, intent(out) :: dsdn(nc) !< J/K/mol - Entropy differential wrpt. mol numbers
    integer, optional, intent(in) :: contribution !< Contribution from ideal (PROP_IDEAL), residual (PROP_RESIDUAL) or both (PROP_OVERALL)
    ! Locals
    real :: s_id, F, F_T, F_V, dsdne(nce), ne(nce)
    real :: sumne, P, dPdV, dPdT, dVdT, zFac, dPdn(nce), dVdn(nce)
    logical :: res, ideal_gas, differentials
    real, pointer :: F_TT_p, F_TV_p, F_VV_p
    real, target :: F_TT, F_TV, F_VV
    real, pointer :: F_Tn_p(:), F_n_p(:), F_Vn_p(:)
    real, target :: F_Tn(nce), F_n(nce), F_Vn(nce)
    real, pointer :: dsdP_id_p, dsdT_id_p, dsdne_id_p(:)
    real, target :: dsdP_id, dsdT_id, dsdne_id(nce)
    !--------------------------------------------------------------------
    res = .true.
    ideal_gas = .true.
    if (present(contribution)) then
      if (contribution == PROP_RESIDUAL) then
        ideal_gas = .false.
      else if (contribution == PROP_IDEAL) then
        res = .false.
      endif
    endif

    F_TV_p => NULL()
    F_VV_p => NULL()
    F_TT_p => NULL()
    F_n_p => NULL()
    F_Tn_p => NULL()
    F_Vn_p => NULL()
    dsdP_id_p => NULL()
    dsdT_id_p => NULL()
    dsdne_id_p => NULL()
    differentials = .false.
    s = 0
    if (present(dsdT) .or. present(dsdp) .or. present(dsdn)) then
      differentials = .true.
      F_TV_p => F_TV
      F_VV_p => F_VV
      if (present(dsdT)) then
        F_TT_p => F_TT
        dsdT = 0
        dsdT_id_p => dsdT_id
      endif
      if (present(dsdP)) then
        dsdP_id_p => dsdP_id
        dsdP = 0
      endif
      if (present(dsdn)) then
        F_n_p => F_n
        F_Tn_p => F_Tn
        F_Vn_p => F_Vn
        dsdne_id_p => dsdne_id
        dsdne = 0
      endif
    endif

    call apparent_to_real_mole_numbers(n,ne)

    if (res) then
      ! Residual contribution
      call Fres_ne(T,v,ne,F=F,F_T=F_T,F_V=F_V,F_n=F_n_p,F_VV=F_VV_p,&
           F_TT=F_TT_p,F_TV=F_TV_p,F_Tn=F_Tn_p,F_Vn=F_Vn_p)

      sumne = sum(ne)
      P = -Rgas*T*F_V + sumne*Rgas*T/V
      zFac = V*P/(sumne*Rgas*T)
      s = Rgas*(-F - T*F_T  + sumne*log(zFac))
      if (differentials) then
        dPdV = -Rgas*T*(F_VV + sumne/V**2)
        dPdT = P/T-Rgas*T*F_TV
        dVdT = -dPdT/dPdV

        if (present(dSdt)) then
          dSdt = dVdt*dPdt - Rgas*(2*F_T + T*F_TT + sumne/T)
        endif

        if (present(dSdp) ) then
          dSdp = sumne*Rgas/P - dVdT
        end if

        if (present(dSdn)) then
          dPdn = Rgas*T*(-F_Vn + 1/V)
          dVdn = -dPdn/dPdV
          dSdne = dVdn*dPdt - Rgas*(F_n + T*F_Tn + 1 - log(zFac))
        endif
      endif
    endif

    if (ideal_gas) then
      dsdt_id = 0
      dsdp_id = 0
      if (.not. res) then
        P = pressure(T,v,n)
      endif
      call idealEntropy_ne(T, P, ne, s_id, dsdT_id_p, dsdP_id_p, dsdne_id_p)
      s = s + s_id
      if (present(dsdt)) then
        dsdt = dsdt + dsdt_id
      endif
      if (present(dsdp)) then
        dsdp = dsdp + dsdp_id
      endif
      if (present(dsdn)) then
        dsdne = dsdne + dsdne_id
      endif
    endif

    if (present(dSdn)) then
      call real_to_apparent_diff(dSdne,dSdn)
    endif
  end subroutine entropy_tvp

end module eosTV

max_stars_repo_name: SINTEF/Thermopack

max_issues_repo_name: SINTEF/Thermopack

max_forks_repo_name: SINTEF/Thermopack


text:
! Please note that more features are available in the Python version

program main
  ! use cmd_line
  use core
  use utils
  implicit none

  character(len=MAX_PATH_LEN) :: config_path, output_path
  real(kind=REAL_KIND), allocatable :: m(:), x(:, :), v(:, :), a(:, :)
  real(kind=REAL_KIND) :: dt, g, min_dist
  integer :: n_steps, n_objs, print_interval, write_interval, ios

  ! Argument processing
  call get_paths(config_path, output_path)

  ! Config processing
  print *, "Using configuration file: ", trim(config_path)
  print *, "Output will be written to: ", trim(output_path)
  call read_config(config_path, x, v, m, dt, n_steps, n_objs, g, min_dist, print_interval, write_interval)
  print *, "n_objs: ", n_objs
  print *, "n_steps: ", n_steps
  print *, "dt: ", dt
  print *, "G: ", g
  print *, "min_dist: ", min_dist
  print *, "print_interval: ", print_interval
  print *, "write_interval: ", write_interval
  call print_arr_2d(x, "x")
  call print_arr_2d(v, "v")
  call print_arr_1d(m, "m")

  allocate(a(DIMS, n_objs))

  ! Output directory processing
  call system("mkdir -p " // trim(output_path), status=ios)
  if (ios /= 0) then
    print *, "Failed to create output directory"
    stop
  end if
  call system("rm -f " // trim(output_path) // "/*.txt", status=ios)
  if (ios /= 0) then
    print *, "Failed to clean output directory"
    stop
  end if

  print *, "Simulating"
  call iterate(x, v, a, m, dt, n_steps, n_objs, g, min_dist, print_interval, write_interval, output_path)
end program main

max_stars_repo_name: AgenttiX/planetary-motion

max_issues_repo_name: AgenttiX/planetary-motion

max_forks_repo_name: AgenttiX/planetary-motion


text:
!*********************************************************************
! This file is part of Polylogarithm.
!
! Polylogarithm is licenced under the MIT License.
!*********************************************************************

!*********************************************************************
!> @brief Clausen function \f$\mathrm{Cl}_2(\theta) = \mathrm{Im}(\mathrm{Li}_2(e^{i\theta}))\f$
!> @param x real angle
!> @return \f$\mathrm{Cl}_2(\theta)\f$
!> @author Alexander Voigt
!> Implemented as economized Padé approximation.
!*********************************************************************

double precision function dcl2(x)
  implicit none
  double precision :: x, y, z, z2, z4, p, q, p0, p1, h, sgn
  double precision, parameter :: PI = 3.14159265358979324D0
  double precision, parameter :: PI2 = 2*PI, PIH = PI/2, PI28 = PI*PI/8
  double precision, parameter :: cp(4) = (/ &
      1.3975782911209635D-2,                &
     -4.4432680257270761D-4,                &
      3.4141174111242951D-6,                &
     -3.7638116201783404D-9                 /)
  double precision, parameter :: cq(4) = (/ &
      1.0000000000000000D+0,                &
     -3.6904397961160525D-2,                &
      3.7342870576106476D-4,                &
     -8.7460760866531179D-7                 /)
  double precision, parameter :: cr(6) = (/ &
      6.4005702446195512D-1,                &
     -2.0641655351338783D-1,                &
      2.4175305223497718D-2,                &
     -1.2355955287855728D-3,                &
      2.5649833551291124D-5,                &
     -1.4783829128773320D-7                 /)
  double precision, parameter :: cs(6) = (/ &
      1.0000000000000000D+0,                &
     -2.5299102015666356D-1,                &
      2.2148751048467057D-2,                &
     -7.8183920462457496D-4,                &
      9.5432542196310670D-6,                &
     -1.8184302880448247D-8                 /)

  sgn = 1

  if (x .lt. 0) then
     x = -x
     sgn = -1
  endif

  if (x .ge. PI2) then
     x = mod(x, PI2)
  endif

  if (x .gt. PI) then
     p0 = 6.28125D0
     p1 = 0.0019353071795864769253D0
     x = (p0 - x) + p1
     sgn = -sgn
  endif

  if (x .eq. 0 .or. x .eq. PI) then
     dcl2 = 0
     return
  endif

  if (x .lt. PIH) then
    y = x*x
    z = y - PI28
    z2 = z*z
    p = cp(1) + z * cp(2) + z2 * (cp(3) + z * cp(4))
    q = cq(1) + z * cq(2) + z2 * (cq(3) + z * cq(4))
    h = x*(1 - log(x) + y*p/q)
  else
    y = PI - x
    z = y*y - PI28
    z2 = z*z
    z4 = z2*z2
    p = cr(1) + z * cr(2) + z2 * (cr(3) + z * cr(4)) + &
        z4 * (cr(5) + z * cr(6))
    q = cs(1) + z * cs(2) + z2 * (cs(3) + z * cs(4)) + &
        z4 * (cs(5) + z * cs(6))
    h = y*p/q
  endif

  dcl2 = sgn*h

end function dcl2

max_stars_repo_name: Expander/polylogarithm

max_issues_repo_name: Expander/polylogarithm

max_forks_repo_name: Expander/polylogarithm


text:
function rnd()
 integer dummmy
 real*4 rnd,ran1
 external ran1
 rnd=ran1(dummy)
end function

function nrnd()
 real,save :: x,r=-1
 real :: nrnd,r1,y
 if(r.ge.0) then
  nrnd=x
  r=-1
 endif
 r=1
 do while(r.ge.1)
  x=2*rnd()-1
  y=2*rnd()-1
  r=sqrt(x*x+y*y)
 enddo
  if(r.ne.0) then
   r1=2*sqrt(-log(r))
   x=r1/r*x
   nrnd=r1/r*y
  else
   nrnd=0.
  endif
end function

function sech2rnd()
 real sech2rnd
 sech2rnd=atanh(1-2*rnd())
end function

subroutine generate_gas(ngas,vdisp,iseed,m,x,y,z,vx,vy,vz,cs)
  use hydrostaticdiskMod
 real, parameter :: pi=3.141592654
 real, allocatable :: rad(:),crad(:),csigtab(:)
 real nrnd,mass,drs,vrot,fipart,rpart,rrpart,zpart,cm,dr
 integer ngas,nrs,bin,bin1,imax,i,j,jm,nfis,nr,ibuf(100)
 real m(ngas),x(ngas),y(ngas),z(ngas),vx(ngas),vy(ngas),vz(ngas),cs(ngas)
 real fr,fz,phi
 character*60 filename
 integer*8 iseed,itmp
 integer*8,parameter :: skip=100

 nrs=1000

 filename='dbh.dat'
 call readharmfile(filename,ibuf)

 call readgasdisk
 nr=nrQ()

 vdisp=0
 vdispt=csQ()
 allocate(rad(0:nr),csigtab(0:nr),crad(0:nrs))
 rad=radQ()
 csigtab=csigtabQ()

 imax=0
 do i=1,nr
  if(csigtab(i).GT.csigtab(imax)) imax=i
!  print*,rad(i),1-csigtab(i)/csigtab(nr)
 enddo
 mass=csigtab(imax)

 dr=rad(imax)/imax
 drs=rad(imax)/nrs
 crad(0)=0
 do i=1,nrs-1
  cm=i*mass/nrs
  do j=0,imax
   if(csigtab(j).lt.cm) jm=j
  enddo
  crad(i)=dr*jm+dr*(cm-csigtab(jm))/(csigtab(jm+1)-csigtab(jm))
 enddo
 crad(nrs)=rad(imax)

 write(0,*) 'cdisp:',vdispt
 write(0,*) 'total mass of gasdisk:',mass,nrs,drs
 write(0,*) 'rmax, mass:',rad(imax),mass/ngas

m=mass/ngas
cs=vdispt
 do i=1,ngas
   itmp=(iseed+i)*skip
   call ran_seed(itmp)
!  bin=(i-1)/nfis+1
!  rpart=sqrt(rnd()*(crad(bin)**2-crad(bin-1)**2)+crad(bin-1)**2)
  rpart=getrad(rnd())
!  fipart=2*pi*(mod(i-1,nfis)+rnd())/nfis
  fipart=2*pi*rnd()
  zpart=getz(rpart,rnd())

  x(i)=rpart*cos(fipart)
  y(i)=rpart*sin(fipart)
  z(i)=zpart
  if(rnd().GT.0.5) z(i)=-z(i)

  call force(rpart,zpart,fr,fz,phi)
  vrot=sqrt(max(0.,-rpart*fr))
!  vrot=sqrt(max(0.,rpart*fr))
! drift correction: checkcheck
!  print*,'zz',rpart,zpart,fr,phi
  vrot=sqrt(max(0.,vrot**2-driftcorfac(rpart)*vdispt**2))
  vx(i)=-vrot*sin(fipart)
  vy(i)=vrot*cos(fipart)
  vz(i)=0

  vx(i)=vx(i)+vdisp*nrnd()
  vy(i)=vy(i)+vdisp*nrnd()
  vz(i)=vz(i)+vdisp*nrnd()

 enddo
end subroutine

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


!#==============================================================================================
!#
!#   Copyright 2020 Abhijit Chatterjee
!#
!#   Licensed under the Apache License, Version 2.0 (the "License");
!#   you may not use this file except in compliance with the License.
!#   You may obtain a copy of the License at
!#
!#       http://www.apache.org/licenses/LICENSE-2.0
!#
!#   Unless required by applicable law or agreed to in writing, software
!#   distributed under the License is distributed on an "AS IS" BASIS,
!#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!#   See the License for the specific language governing permissions and
!#   limitations under the License.
!#==============================================================================================

!#```````````````````````````````````````````````````````````````````````````````````````````````````````````````````````````````````````
!#Script name: tpmdvariables.f90
!#
!#Purpose: Basic TPMD variables
!#
!#Authors: Abhijit Chatterjee (abhijit@che.iitb.ac.in)
!#
!#Date of Modification: 15 May 2020
!#.......................................................................................................................................

MODULE TPMDVariables
   IMPLICIT NONE

   DOUBLE PRECISION :: kB=8.617343000e-5
   TYPE TPMDData
      INTEGER :: npts=0
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: time
      INTEGER, DIMENSION(:), ALLOCATABLE :: fstate
      DOUBLE PRECISION :: tau=0.,T0=0.,deltaT=0.,Tmax=0.
   END TYPE TPMDData
   TYPE TPMDStagePartition
      INTEGER :: Nmax=0,m=0
      DOUBLE PRECISION :: invTna_sum=0.d0
      INTEGER, DIMENSION(:,:), ALLOCATABLE :: StagewiseBreakupState
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: StagewiseBreakupTime
      INTEGER, DIMENSION(:), ALLOCATABLE :: count,StagewisePoints
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: theta
   END TYPE TPMDStagePartition

   CONTAINS

   SUBROUTINE InitializeRandomSeed(iseed)
      IMPLICIT NONE
      INTEGER :: iseed,n
      integer, allocatable :: seed(:)

      call random_seed(size = n)
      allocate(seed(n))
      seed=iseed
      call random_seed(put=seed)
   END SUBROUTINE InitializeRandomSeed
   !xxxxxxxxxxxxx
   SUBROUTINE CreateTPMDDataset(d,tau,T0,deltaT,Tmax)
      IMPLICIT NONE
      TYPE(TPMDData), POINTER :: d
      DOUBLE PRECISION :: tau,T0,deltaT,Tmax

      ALLOCATE(d)
      d%tau=tau
      d%T0=T0
      d%deltaT=deltaT
      d%Tmax=Tmax
   END SUBROUTINE CreateTPMDDataset
   !xxxxxxxxxxxxx
   SUBROUTINE DeleteTPMDDataset(d)
      IMPLICIT NONE
      TYPE(TPMDData), POINTER :: d

      IF (ALLOCATED(d%time)) DEALLOCATE(d%time)
      IF (ALLOCATED(d%fstate)) DEALLOCATE(d%fstate)
      IF (ASSOCIATED(d)) DEALLOCATE(d)
   END SUBROUTINE DeleteTPMDDataset
END MODULE TPMDVariables
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx