code: # Support for python2
  from __future__ import print_function
  
  import sys
  from ctypes import *
  from os import path
  import numpy as np
  from . import plotutils, utils, platform_specifics
  
  if utils.gcc_major_version_greater_than(7):
      c_len_type = c_size_t # c_size_t on GCC > 7
  else:
      c_len_type = c_int
  
  
  class thermopack(object):
      """
      Interface to thermopack
      """
      def __init__(self):
          """
          Load libthermopack.(so/dll) and initialize function pointers
          """
          pf_specifics = platform_specifics.get_platform_specifics()
          self.prefix = pf_specifics["prefix"]
          self.module = pf_specifics["module"]
          self.postfix = pf_specifics["postfix"]
          self.postfix_nm = pf_specifics["postfix_no_module"]
          dyn_lib_path = path.join(path.dirname(__file__), pf_specifics["dyn_lib"])
          self.tp = cdll.LoadLibrary(dyn_lib_path)
  
          # Set phase flags
          self.s_get_phase_flags = self.tp.get_phase_flags_c
          self.get_phase_flags()
  
          # Model control
          self.s_add_eos = getattr(self.tp, self.get_export_name("thermopack_var", "add_eos"))
          self.s_delete_eos = getattr(self.tp, self.get_export_name("thermopack_var", "delete_eos"))
          self.s_activate_model = getattr(self.tp, self.get_export_name("thermopack_var", "activate_model"))
  
          # Information
          self.s_get_model_id = getattr(self.tp, self.get_export_name("thermopack_var", "get_eos_identification"))
  
          # Init methods
          self.eoslibinit_init_thermo = getattr(self.tp, self.get_export_name("eoslibinit", "init_thermo"))
          self.Rgas = c_double.in_dll(self.tp, self.get_export_name("thermopack_constants", "rgas")).value
          self.nc = None
          self.minimum_temperature_c = c_double.in_dll(self.tp, self.get_export_name("thermopack_constants", "tptmin"))
          self.minimum_pressure_c = c_double.in_dll(self.tp, self.get_export_name("thermopack_constants", "tppmin"))
          self.solideos_solid_init = getattr(self.tp, self.get_export_name("solideos", "solid_init"))
          self.eoslibinit_init_volume_translation = getattr(self.tp, self.get_export_name("eoslibinit", "init_volume_translation"))
          self.eoslibinit_redefine_critical_parameters = getattr(self.tp, self.get_export_name("eoslibinit", "redefine_critical_parameters"))
  
          # Eos interface
          self.s_eos_specificvolume = getattr(self.tp, self.get_export_name("eos", "specificvolume"))
          self.s_eos_zfac = getattr(self.tp, self.get_export_name("eos", "zfac"))
          self.s_eos_thermo = getattr(self.tp, self.get_export_name("eos", "thermo"))
          self.s_eos_entropy = getattr(self.tp, self.get_export_name("eos", "entropy"))
          self.s_eos_enthalpy = getattr(self.tp, self.get_export_name("eos", "enthalpy"))
          self.s_eos_compmoleweight = getattr(self.tp, self.get_export_name("eos", "compmoleweight"))
          self.s_eos_idealenthalpysingle = getattr(self.tp, self.get_export_name("eos", "idealenthalpysingle"))
  
          # Speed of sound
          #self.sos_singlePhaseSpeedOfSound = getattr(self.tp, '__speed_of_sound_MOD_singlephasespeedofsound')
          self.s_sos_sound_velocity_2ph = getattr(self.tp, self.get_export_name("speed_of_sound", "sound_velocity_2ph"))
  
          # Component info
          self.s_compdata_compindex = getattr(self.tp, self.get_export_name("compdata", "comp_index_active"))
          self.s_compdata_compname = getattr(self.tp, self.get_export_name("compdata", "comp_name_active"))
  
          # Flashes
          self.s_set_ph_tolerance = getattr(self.tp, self.get_export_name("ph_solver", "setphtolerance"))
          self.s_twophasetpflash = getattr(self.tp, self.get_export_name("tp_solver", "twophasetpflash"))
          self.s_psflash_twophase = getattr(self.tp, self.get_export_name("ps_solver", "twophasepsflash"))
          #self.tpflash_multiphase = getattr(self.tp, '__mp_tp_solver_MOD_mp_flash_tp')
          self.s_uvflash_twophase = getattr(self.tp, self.get_export_name("uv_solver", "twophaseuvflash"))
          self.s_phflash_twophase = getattr(self.tp, self.get_export_name("ph_solver", "twophasephflash"))
          #self.s_svflash_twophase = getattr(self.tp, self.get_export_name("sv_solver", "twophasesvflash"))
          self.s_guess_phase = getattr(self.tp, self.get_export_name("thermo_utils", "guessphase"))
  
          # TV interfaces
          self.s_internal_energy_tv = getattr(self.tp, self.get_export_name("eostv", "internal_energy"))
          self.s_entropy_tv = getattr(self.tp, self.get_export_name("eostv", "entropytv"))
          self.s_pressure_tv = getattr(self.tp, self.get_export_name("eostv", "pressure"))
          self.s_lnphi_tv = getattr(self.tp, self.get_export_name("eostv", "thermotv"))
          self.s_enthalpy_tv = getattr(self.tp, self.get_export_name("eostv", "enthalpytv"))
          self.s_helmholtz_energy = getattr(self.tp, self.get_export_name("eostv", "free_energy"))
          self.s_chempot = getattr(self.tp, self.get_export_name("eostv", "chemical_potential"))
  
          # Saturation properties
          self.s_bubble_t = getattr(self.tp, self.get_export_name("saturation", "safe_bubt"))
          self.s_bubble_p = getattr(self.tp, self.get_export_name("saturation", "safe_bubp"))
          self.s_dew_t = getattr(self.tp, self.get_export_name("saturation", "safe_dewt"))
          self.s_dew_p = getattr(self.tp, self.get_export_name("saturation", "safe_dewp"))
          self.s_envelope_plot = getattr(self.tp, self.get_export_name("saturation_curve", "envelopeplot"))
          self.s_binary_plot = getattr(self.tp, self.get_export_name("binaryplot", "vllebinaryxy"))
          self.s_global_binary_plot = getattr(self.tp, self.get_export_name("binaryplot", "global_binary_plot"))
          self.s_get_bp_term = getattr(self.tp, self.get_export_name("binaryplot", "get_bp_term"))
          self.s_solid_envelope_plot = getattr(self.tp, self.get_export_name("solid_saturation", "solidenvelopeplot"))
          self.s_isotherm = getattr(self.tp, self.get_export_name("isolines", "isotherm"))
          self.s_isobar = getattr(self.tp, self.get_export_name("isolines", "isobar"))
          self.s_isenthalp = getattr(self.tp, self.get_export_name("isolines", "isenthalp"))
          self.s_isentrope = getattr(self.tp, self.get_export_name("isolines", "isentrope"))
          # Stability
          self.s_crit_tv = getattr(self.tp, self.get_export_name("critical", "calccriticaltv"))
  
          # Virials
          self.s_virial_coeffcients = getattr(self.tp, self.get_export_name("eostv", "virial_coefficients"))
          self.s_second_virial_matrix = getattr(self.tp, self.get_export_name("eostv", "secondvirialcoeffmatrix"))
          self.s_binary_third_virial_matrix = getattr(self.tp, self.get_export_name("eostv", "binarythirdvirialcoeffmatrix"))
  
          self.add_eos()
  
      def __del__(self):
          """Delete FORTRAN memory allocated for this instance"""
          self.delete_eos()
  
      def activate(self):
          """Activate this instance of thermopack parameters for calculation"""
          self.s_activate_model.argtypes = [POINTER( c_int )]
          self.s_activate_model.restype = None
          self.s_activate_model(self.model_index_c)
  
      def add_eos(self):
          """Allocate FORTRAN memory for this class instance"""
          self.s_add_eos.argtypes = None
          self.s_add_eos.restype = c_int
          self.model_index_c = c_int(self.s_add_eos())
  
      def delete_eos(self):
          """de-allocate FORTRAN memory for this class instance"""
          self.activate()
          self.s_delete_eos.argtypes = [POINTER( c_int )]
          self.s_delete_eos.restype = None
          self.s_delete_eos(self.model_index_c)
  
      def get_model_id(self):
          """Get model identification
  
          Returns:
              str: Eos name
          """
          self.activate()
  
          eosid_len = 40
          eosid_c = c_char_p(b" " * eosid_len)
          eosid_len_c = c_len_type(eosid_len)
          self.s_get_model_id.argtypes = [c_char_p, c_len_type]
          self.s_get_model_id.restype = None
          self.s_get_model_id(eosid_c, eosid_len_c)
  
          eosid = eosid_c.value.decode('ascii').strip()
          return eosid
  
      def get_export_name(self, module, method):
          """Generate library export name based on module and method name
  
          Args:
              module (str): Name of module
              method (str): Name of method
  
          Returns:
              str: Library export name
          """
          if len(module) > 0:
              export_name = self.prefix + module + self.module + method + self.postfix
          else:
              export_name = method + self.postfix_nm
          return export_name
  
      #################################
      # Init
      #################################
  
      def init_thermo(self, eos, mixing, alpha, comps, nphases,
                      liq_vap_discr_method=None, csp_eos=None, csp_ref_comp=None,
                      kij_ref="Default", alpha_ref="Default", saft_ref="Default",
                      b_exponent=None, TrendEosForCp=None, cptype=None,
                      silent=None):
          """Initialize thermopack
  
          Args:
              eos (str): Equation of state
              mixing (str): Mixture model for cubic eos
              alpha (str): Alpha formulations for cubic EOS
              comps (string): Comma separated list of components
              nphases (int): Maximum number of phases considered during multi-phase flash calculations
              liq_vap_discr_method (int, optional): Method to discriminate between liquid and vapor in case of an undefined single phase. Defaults to None.
              csp_eos (str, optional): Corrensponding state equation. Defaults to None.
              csp_ref_comp (str, optional): CSP reference component. Defaults to None.
              kij_ref (str, optional): Data set identifiers. Defaults to "Default".
              alpha_ref (str, optional): Data set identifiers. Defaults to "Default".
              saft_ref (str, optional): Data set identifiers. Defaults to "Default".
              b_exponent (float, optional): Exponent used in co-volume mixing. Defaults to None.
              TrendEosForCp (str, optional): Option to init trend for ideal gas properties. Defaults to None.
              cptype (int array, optional): Equation type number for Cp. Defaults to None.
              silent (bool, optional): Supress messages during init?. Defaults to None.
          """
          self.activate()
          self.nc = max(len(comps.split(" ")), len(comps.split(",")))
  
          null_pointer = POINTER(c_int)()
          eos_c = c_char_p(eos.encode('ascii'))
          eos_len = c_len_type(len(eos))
          mixing_c = c_char_p(mixing.encode('ascii'))
          mixing_len = c_len_type(len(mixing))
          alpha_c = c_char_p(alpha.encode('ascii'))
          alpha_len = c_len_type(len(alpha))
          comp_string_c = c_char_p(comps.encode('ascii'))
          comp_string_len = c_len_type(len(comps))
          nphases_c = c_int(nphases)
          if liq_vap_discr_method is None:
              liq_vap_discr_method_c = null_pointer
          else:
              liq_vap_discr_method_c = POINTER(c_int)(c_int(liq_vap_discr_method))
          if csp_eos is None:
              csp_eos_c = c_char_p()
              csp_eos_len = c_len_type(0)
          else:
              csp_eos_c = c_char_p(csp_eos.encode('ascii'))
              csp_eos_len = c_len_type(len(csp_eos))
          if csp_ref_comp is None:
              csp_ref_comp_c = c_char_p()
              csp_ref_comp_len = c_len_type(0)
          else:
              csp_ref_comp_c = c_char_p(csp_ref_comp.encode('ascii'))
              csp_ref_comp_len = c_len_type(len(csp_ref_comp))
          kij_ref_len = c_len_type(len(kij_ref))
          kij_ref_c = c_char_p(kij_ref.encode('ascii'))
          alpha_ref_len = c_len_type(len(alpha_ref))
          alpha_ref_c = c_char_p(alpha_ref.encode('ascii'))
          saft_ref_len = c_len_type(len(saft_ref))
          saft_ref_c = c_char_p(saft_ref.encode('ascii'))
          if b_exponent is None:
              b_exponent_c = POINTER(c_double)()
          else:
              b_exponent_c = POINTER(c_double)(c_double(b_exponent))
          if TrendEosForCp is None:
              TrendEosForCp_c = c_char_p()
              TrendEosForCp_len = c_len_type(0)
          else:
              TrendEosForCp_c = c_char_p(TrendEosForCp.encode('ascii'))
              TrendEosForCp_len = c_len_type(len(TrendEosForCp))
          if cptype is None:
              cptype_c = null_pointer
          else:
              cptype_c = (c_int * self.nc)(*cptype)
  
          if silent is None:
              silent_c = null_pointer
          else:
              if silent:
                  silent_int = 1
              else:
                  silent_int = 0
              silent_c = POINTER(c_int)(c_int(silent_int))
  
          self.eoslibinit_init_thermo.argtypes = [c_char_p,
                                                  c_char_p,
                                                  c_char_p,
                                                  c_char_p,
                                                  POINTER( c_int ),
                                                  POINTER( c_int ),
                                                  c_char_p,
                                                  c_char_p,
                                                  c_char_p,
                                                  c_char_p,
                                                  c_char_p,
                                                  POINTER( c_double ),
                                                  c_char_p,
                                                  POINTER( c_int ),
                                                  POINTER( c_int ),
                                                  c_len_type, c_len_type,
                                                  c_len_type, c_len_type,
                                                  c_len_type, c_len_type,
                                                  c_len_type, c_len_type,
                                                  c_len_type, c_len_type]
  
  
          self.eoslibinit_init_thermo.restype = None
  
          self.eoslibinit_init_thermo(eos_c,
                                      mixing_c,
                                      alpha_c,
                                      comp_string_c,
                                      byref(nphases_c),
                                      liq_vap_discr_method_c,
                                      csp_eos_c,
                                      csp_ref_comp_c,
                                      kij_ref_c,
                                      alpha_ref_c,
                                      saft_ref_c,
                                      b_exponent_c,
                                      TrendEosForCp_c,
                                      cptype_c,
                                      silent_c,
                                      eos_len,
                                      mixing_len,
                                      alpha_len,
                                      comp_string_len,
                                      csp_eos_len,
                                      csp_ref_comp_len,
                                      kij_ref_len,
                                      alpha_ref_len,
                                      saft_ref_len,
                                      TrendEosForCp_len)
  
      def init_peneloux_volume_translation(self, parameter_reference="Default"):
          """Initialialize Peneloux volume translations
  
          Args:
              parameter_reference (str): String defining parameter set, Defaults to "Default"
          """
          self.activate()
          volume_trans_model = "PENELOUX"
          volume_trans_model_c = c_char_p(volume_trans_model.encode('ascii'))
          volume_trans_model_len = c_len_type(len(volume_trans_model))
          ref_string_c = c_char_p(parameter_reference.encode('ascii'))
          ref_string_len = c_len_type(len(parameter_reference))
          self.eoslibinit_init_volume_translation.argtypes = [c_char_p,
                                                              c_char_p,
                                                              c_len_type,
                                                              c_len_type]
  
          self.eoslibinit_init_volume_translation.restype = None
  
          self.eoslibinit_init_volume_translation(volume_trans_model_c,
                                                  ref_string_c,
                                                  volume_trans_model_len,
                                                  ref_string_len)
  
      def redefine_critical_parameters(self, silent=True):
          """Recalculate critical properties of pure fluids
  
          Args:
              silent (bool): Ignore warnings? Defaults to True
          """
          self.activate()
          if silent:
              silent_c = c_int(1)
          else:
              silent_c = c_int(0)
  
          self.eoslibinit_redefine_critical_parameters.argtypes = [ POINTER( c_int ) ]
  
          self.eoslibinit_redefine_critical_parameters.restype = None
  
          self.eoslibinit_redefine_critical_parameters(byref(silent_c))
  
  
      #################################
      # Solids
      #################################
  
      def init_solid(self, scomp):
          """Initialize pure solid
  
          Args:
              scomp (str): Component name
          """
          self.activate()
          scomp_c = c_char_p(scomp.encode('ascii'))
          scomp_len = c_len_type(len(scomp))
          self.solideos_solid_init.argtypes = [c_char_p, c_len_type]
          self.solideos_solid_init.restype = None
          self.solideos_solid_init(scomp_c, scomp_len)
  
      #################################
      # Utility
      #################################
  
      def getcompindex(self, comp):
          """Get component index
  
          Args:
              comp (str): Component name
  
          Returns:
              int: Component FORTRAN index
          """
          self.activate()
          comp_c = c_char_p(comp.encode('ascii'))
          comp_len = c_len_type(len(comp))
          self.s_compdata_compindex.argtypes = [c_char_p, c_len_type]
          self.s_compdata_compindex.restype = c_int
          idx = self.s_compdata_compindex(comp_c, comp_len)
          return idx
  
      def get_comp_name(self, index):
          """Get component name
  
          Args:
              int: Component FORTRAN index
  
          Returns:
              comp (str): Component name
          """
          self.activate()
          comp_len = 40
          comp_c = c_char_p(b" " * comp_len)
          comp_len_c = c_len_type(comp_len)
          index_c = c_int(index)
          self.s_compdata_compname.argtypes = [POINTER(c_int), c_char_p, c_len_type]
          self.s_compdata_compname.restype = None
          self.s_compdata_compname(byref(index_c), comp_c, comp_len_c)
          compname = comp_c.value.decode('ascii').strip()
          return compname
  
      def compmoleweight(self, comp):
          """Get component mole weight (g/mol)
  
          Args:
              comp (int): Component FORTRAN index
  
          Returns:
              float: Component mole weight (g/mol)
          """
          self.activate()
          comp_c = c_int(comp)
          self.s_eos_compmoleweight.argtypes = [POINTER( c_int )]
          self.s_eos_compmoleweight.restype = c_double
          mw_i = self.s_eos_compmoleweight(byref(comp_c))
          return mw_i
  
      def get_phase_flags(self):
          """Get phase identifiers used by thermopack
  
          Returns:
              int: Phase int  identifiers
          """
          iTWOPH = c_int()
          iLIQPH = c_int()
          iVAPPH = c_int()
          iMINGIBBSPH = c_int()
          iSINGLEPH = c_int()
          iSOLIDPH = c_int()
          iFAKEPH = c_int()
  
          self.s_get_phase_flags.argtypes = [POINTER( c_int ),
                                             POINTER( c_int ),
                                             POINTER( c_int ),
                                             POINTER( c_int ),
                                             POINTER( c_int ),
                                             POINTER( c_int ),
                                             POINTER( c_int )]
          self.s_get_phase_flags.restype = None
          self.s_get_phase_flags(byref(iTWOPH),
                                 byref(iLIQPH),
                                 byref(iVAPPH),
                                 byref(iMINGIBBSPH),
                                 byref(iSINGLEPH),
                                 byref(iSOLIDPH),
                                 byref(iFAKEPH))
          self.TWOPH = iTWOPH.value
          self.LIQPH = iLIQPH.value
          self.VAPPH = iVAPPH.value
          self.MINGIBBSPH = iMINGIBBSPH.value
          self.SINGLEPH = iSINGLEPH.value
          self.SOLIDPH = iSOLIDPH.value
          self.FAKEPH = iFAKEPH.value
  
      def get_phase_type(self, i_phase):
          """Get phase type
  
          Args:
              i_phase (int): Phase flag returned by thermopack
  
          Returns:
              str: Phase type
          """
          phase_string_list = ["TWO_PHASE", "LIQUID", "VAPOR", "MINIMUM_GIBBS", "SINGLE", "SOLID", "FAKE"]
          return phase_string_list[i_phase]
  
      def set_tmin(self, temp):
          """Set minimum temperature in Thermopack. Used to limit search
          domain for numerical solvers.
  
          Args:
              temp (float): Temperature (K)
          """
          self.minimum_temperature_c.value = temp
  
      def get_tmin(self):
          """Get minimum temperature in Thermopack. Used to limit search
          domain for numerical solvers.
  
          Returns:
              float: Temperature (K)
          """
          temp = self.minimum_temperature_c.value
          return temp
  
      def set_pmin(self, press):
          """Get minimum pressure in Thermopack. Used to limit search
          domain for numerical solvers.
  
          Args:
              press (float): Pressure (Pa)
          """
          self.minimum_pressure_c.value = press
  
      #################################
      # Phase properties
      #################################
  
      def specific_volume(self, temp, press, x, phase, dvdt=None, dvdp=None, dvdn=None):
          """ Calculate single-phase specific volume
              Note that the order of the output match the default order of input for the differentials.
              Note further that dvdt, dvdp and dvdn only are flags to enable calculation.
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (array_like): Molar composition
              phase (int): Calcualte root for specified phase
              dvdt (logical, optional): Calculate volume differentials with respect to temperature while pressure and composition are held constant. Defaults to None.
              dvdp (logical, optional): Calculate volume differentials with respect to pressure while temperature and composition are held constant. Defaults to None.
              dvdn (logical, optional): Calculate volume differentials with respect to mol numbers while pressure and temperature are held constant. Defaults to None.
  
          Returns:
              float: Specific volume (m3/mol), and optionally differentials
          """
          self.activate()
          null_pointer = POINTER(c_double)()
  
          temp_c = c_double(temp)
          press_c = c_double(press)
          x_c = (c_double * len(x))(*x)
          phase_c = c_int(phase)
          v_c = c_double(0.0)
  
          if dvdt is None:
              dvdt_c = null_pointer
          else:
              dvdt_c = POINTER(c_double)(c_double(0.0))
          if dvdp is None:
              dvdp_c = null_pointer
          else:
              dvdp_c = POINTER(c_double)(c_double(0.0))
          if dvdn is None:
              dvdn_c = null_pointer
          else:
              dvdn_c = (c_double * len(x))(0.0)
  
          self.s_eos_specificvolume.argtypes = [POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_int ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double )]
  
          self.s_eos_specificvolume.restype = None
  
          self.s_eos_specificvolume(byref(temp_c),
                                    byref(press_c),
                                    x_c,
                                    byref(phase_c),
                                    byref(v_c),
                                    dvdt_c,
                                    dvdp_c,
                                    dvdn_c)
          return_tuple = (v_c.value, )
          if not dvdt is None:
              return_tuple += (dvdt_c[0], )
          if not dvdp is None:
              return_tuple += (dvdp_c[0], )
          if not dvdn is None:
              return_tuple += (np.array(dvdn_c), )
  
          return return_tuple
  
      def zfac(self,temp,press,x,phase,dzdt=None,dzdp=None,dzdn=None):
          """ Calculate single-phase compressibility
              Note that the order of the output match the default order of input for the differentials.
              Note further that dzdt, dzdp and dzdn only are flags to enable calculation.
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (array_like): Molar composition
              phase (int): Calcualte root for specified phase
              dzdt (logical, optional): Calculate compressibility differentials with respect to temperature while pressure and composition are held constant. Defaults to None.
              dzdp (logical, optional): Calculate compressibility differentials with respect to pressure while temperature and composition are held constant. Defaults to None.
              dzdn (logical, optional): Calculate compressibility differentials with respect to mol numbers while pressure and temperature are held constant. Defaults to None.
  
          Returns:
              float: Compressibility (-), and optionally differentials
          """
          self.activate()
          null_pointer = POINTER(c_double)()
  
          temp_c = c_double(temp)
          press_c = c_double(press)
          x_c = (c_double * len(x))(*x)
          phase_c = c_int(phase)
          z_c = c_double(0.0)
  
          if dzdt is None:
              dzdt_c = null_pointer
          else:
              dzdt_c = POINTER(c_double)(c_double(0.0))
          if dzdp is None:
              dzdp_c = null_pointer
          else:
              dzdp_c = POINTER(c_double)(c_double(0.0))
          if dzdn is None:
              dzdn_c = null_pointer
          else:
              dzdn_c = (c_double * len(x))(0.0)
  
          self.s_eos_zfac.argtypes = [POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_int ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double )]
  
          self.s_eos_zfac.restype = None
  
          self.s_eos_zfac(byref(temp_c),
                          byref(press_c),
                          x_c,
                          byref(phase_c),
                          byref(z_c),
                          dzdt_c,
                          dzdp_c,
                          dzdn_c)
          return_tuple = (z_c.value, )
          if not dzdt is None:
              return_tuple += (dzdt_c[0], )
          if not dzdp is None:
              return_tuple += (dzdp_c[0], )
          if not dzdn is None:
              return_tuple += (np.array(dzdn_c), )
  
          return return_tuple
  
      def thermo(self,temp,press,x,phase,dlnfugdt=None,dlnfugdp=None,
                 dlnfugdn=None,ophase=None,v=None):
          """ Calculate logarithm of fugacity coefficient given composition,
          temperature and pressure.
          Note that the order of the output match the default order of input for the differentials.
          Note further that dlnfugdt, dlnfugdp, dlnfugdn and ophase only are flags to enable calculation.
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (array_like): Molar composition (.)
              phase (int): Calcualte root for specified phase
              dlnfugdt (logical, optional): Calculate fugacity coefficient differentials with respect to temperature while pressure and composition are held constant. Defaults to None.
              dlnfugdp (logical, optional): Calculate fugacity coefficient differentials with respect to pressure while temperature and composition are held constant. Defaults to None.
              dlnfugdn (logical, optional): Calculate fugacity coefficient differentials with respect to mol numbers while pressure and temperature are held constant. Defaults to None.
              ophase (int, optional): Phase flag. Only set when phase=MINGIBBSPH.
              v (float, optional): Specific volume (m3/mol)
          Returns:
              ndarray: fugacity coefficient (-), and optionally differentials
          """
          self.activate()
          null_pointer = POINTER(c_double)()
          temp_c = c_double(temp)
          press_c = c_double(press)
          x_c = (c_double * len(x))(*x)
          phase_c = c_int(phase)
          lnfug_c = (c_double * len(x))(0.0)
  
          if dlnfugdt is None:
              dlnfugdt_c = null_pointer
          else:
              dlnfugdt_c = (c_double * len(x))(0.0)
          if dlnfugdp is None:
              dlnfugdp_c = null_pointer
          else:
              dlnfugdp_c = (c_double * len(x))(0.0)
          if dlnfugdn is None:
              dlnfugdn_c = null_pointer
          else:
              dlnfugdn_c = (c_double * len(x)**2)(0.0)
          if ophase is None:
              ophase_c = POINTER(c_int)()
          else:
              ophase_c = POINTER(c_int)(c_int(0))
          metaExtremum_c = POINTER(c_int)()
  
          if v is None:
              v_c = null_pointer
          else:
              v_c = POINTER(c_double)(c_double(0.0))
  
          self.s_eos_thermo.argtypes = [POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_int ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_int ),
                                        POINTER( c_int ),
                                        POINTER( c_double )]
  
          self.s_eos_thermo.restype = None
  
          self.s_eos_thermo(byref(temp_c),
                            byref(press_c),
                            x_c,
                            byref(phase_c),
                            lnfug_c,
                            dlnfugdt_c,
                            dlnfugdp_c,
                            dlnfugdn_c,
                            ophase_c,
                            metaExtremum_c,
                            v_c)
  
          return_tuple = (np.array(lnfug_c), )
          if not dlnfugdt is None:
              return_tuple += (np.array(dlnfugdt_c), )
          if not dlnfugdp is None:
              return_tuple += (np.array(dlnfugdp_c), )
          if not dlnfugdn is None:
              dlnfugdn_r = np.zeros((len(x),len(x)))
              for i in range(len(x)):
                  for j in range(len(x)):
                      dlnfugdn_r[i][j] = dlnfugdn_c[i+j*len(x)]
              return_tuple += (dlnfugdn_r, )
          if not ophase is None:
              return_tuple += (ophase_c[0], )
          if not v is None:
              return_tuple += (v_c[0], )
  
          return return_tuple
  
      def enthalpy(self,temp,press,x,phase,dhdt=None,dhdp=None,dhdn=None):
          """ Calculate specific single-phase enthalpy
              Note that the order of the output match the default order of input for the differentials.
              Note further that dhdt, dhdp and dhdn only are flags to enable calculation.
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (array_like): Molar composition
              phase (int): Calcualte root for specified phase
              dhdt (logical, optional): Calculate enthalpy differentials with respect to temperature while pressure and composition are held constant. Defaults to None.
              dhdp (logical, optional): Calculate enthalpy differentials with respect to pressure while temperature and composition are held constant. Defaults to None.
              dhdn (logical, optional): Calculate enthalpy differentials with respect to mol numbers while pressure and temperature are held constant. Defaults to None.
  
          Returns:
              float: Specific enthalpy (J/mol), and optionally differentials
          """
          self.activate()
          null_pointer = POINTER(c_double)()
  
          temp_c = c_double(temp)
          press_c = c_double(press)
          x_c = (c_double * len(x))(*x)
          phase_c = c_int(phase)
          h_c = c_double(0.0)
  
          if dhdt is None:
              dhdt_c = null_pointer
          else:
              dhdt_c = POINTER(c_double)(c_double(0.0))
          if dhdp is None:
              dhdp_c = null_pointer
          else:
              dhdp_c = POINTER(c_double)(c_double(0.0))
          if dhdn is None:
              dhdn_c = null_pointer
          else:
              dhdn_c = (c_double * len(x))(0.0)
  
          residual_c = POINTER(c_int)()
          self.s_eos_enthalpy.argtypes = [POINTER( c_double ),
                                          POINTER( c_double ),
                                          POINTER( c_double ),
                                          POINTER( c_int ),
                                          POINTER( c_double ),
                                          POINTER( c_double ),
                                          POINTER( c_double ),
                                          POINTER( c_double ),
                                          POINTER( c_int )]
  
          self.s_eos_enthalpy.restype = None
  
          self.s_eos_enthalpy(byref(temp_c),
                              byref(press_c),
                              x_c,
                              byref(phase_c),
                              byref(h_c),
                              dhdt_c,
                              dhdp_c,
                              dhdn_c,
                              residual_c)
  
          return_tuple = (h_c.value, )
          if not dhdt is None:
              return_tuple += (dhdt_c[0], )
          if not dhdp is None:
              return_tuple += (dhdp_c[0], )
          if not dhdn is None:
              return_tuple += (np.array(dhdn_c), )
  
          return return_tuple
  
      def entropy(self,temp,press,x,phase,dsdt=None,dsdp=None,dsdn=None):
          """ Calculate specific single-phase entropy
              Note that the order of the output match the default order of input for the differentials.
              Note further that dsdt, dhsp and dsdn only are flags to enable calculation.
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (array_like): Molar composition
              phase (int): Calcualte root for specified phase
              dsdt (logical, optional): Calculate entropy differentials with respect to temperature while pressure and composition are held constant. Defaults to None.
              dsdp (logical, optional): Calculate entropy differentials with respect to pressure while temperature and composition are held constant. Defaults to None.
              dsdn (logical, optional): Calculate entropy differentials with respect to mol numbers while pressure and temperature are held constant. Defaults to None.
  
          Returns:
              float: Specific entropy (J/mol/K), and optionally differentials
          """
          self.activate()
          null_pointer = POINTER(c_double)()
  
          temp_c = c_double(temp)
          press_c = c_double(press)
          x_c = (c_double * len(x))(*x)
          phase_c = c_int(phase)
          s_c = c_double(0.0)
  
          if dsdt is None:
              dsdt_c = null_pointer
          else:
              dsdt_c = POINTER(c_double)(c_double(0.0))
          if dsdp is None:
              dsdp_c = null_pointer
          else:
              dsdp_c = POINTER(c_double)(c_double(0.0))
          if dsdn is None:
              dsdn_c = null_pointer
          else:
              dsdn_c = (c_double * len(x))(0.0)
          residual_c = POINTER(c_int)()
  
          self.s_eos_entropy.argtypes = [POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_int ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_int )]
  
          self.s_eos_entropy.restype = None
  
          self.s_eos_entropy(byref(temp_c),
                             byref(press_c),
                             x_c,
                             byref(phase_c),
                             byref(s_c),
                             dsdt_c,
                             dsdp_c,
                             dsdn_c,
                             residual_c)
          return_tuple = (s_c.value, )
          if not dsdt is None:
              return_tuple += (dsdt_c[0], )
          if not dsdp is None:
              return_tuple += (dsdp_c[0], )
          if not dsdn is None:
              return_tuple += (np.array(dsdn_c), )
  
          return return_tuple
  
      def idealenthalpysingle(self,temp,press,j,dhdt=None,dhdp=None):
          """ Calculate specific ideal enthalpy
              Note that the order of the output match the default order of input for the differentials.
              Note further that dhdt, and dhdp only are flags to enable calculation.
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (array_like): Molar composition
              phase (int): Calcualte root for specified phase
              dhdt (logical, optional): Calculate ideal enthalpy differentials with respect to temperature while pressure and composition are held constant. Defaults to None.
              dhdp (logical, optional): Calculate ideal enthalpy differentials with respect to pressure while temperature and composition are held constant. Defaults to None.
  
          Returns:
              float: Specific ideal enthalpy (J/mol), and optionally differentials
          """
          self.activate()
          null_pointer = POINTER(c_double)()
  
          temp_c = c_double(temp)
          press_c = c_double(press)
          j_c = c_int(j)
          h_c = c_double(0.0)
  
          if dhdt is None:
              dhdt_c = null_pointer
          else:
              dhdt_c = POINTER(c_double)(c_double(0.0))
          if dhdp is None:
              dhdp_c = null_pointer
          else:
              dhdp_c = POINTER(c_double)(c_double(0.0))
  
          self.s_eos_idealenthalpysingle.argtypes = [POINTER( c_double ),
                                                     POINTER( c_double ),
                                                     POINTER( c_int ),
                                                     POINTER( c_double ),
                                                     POINTER( c_double ),
                                                     POINTER( c_double )]
  
          self.s_eos_idealenthalpysingle.restype = None
  
          self.s_eos_idealenthalpysingle(byref(temp_c),
                                         byref(press_c),
                                         byref(j_c),
                                         byref(h_c),
                                         dhdt_c,
                                         dhdp_c)
          return_tuple = (h_c.value, )
          if not dhdt is None:
              return_tuple += (dhdt_c[0], )
          if not dhdp is None:
              return_tuple += (dhdp_c[0], )
  
          return return_tuple
  
      def speed_of_sound(self,temp,press,x,y,z,betaV,betaL,phase):
          """Calculate speed of sound for single phase or two phase mixture assuming
          mechanical, thermal and chemical equilibrium.
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (array_like): Liquid molar composition
              y (array_like): Gas molar composition
              z (array_like): Overall molar composition
              betaV (float): Molar gas phase fraction
              betaL (float): Molar liquid phase fraction
              phase (int): Calcualte root for specified phase
  
          Returns:
              float: Speed of sound (m/s)
          """
          self.activate()
          temp_c = c_double(temp)
          press_c = c_double(press)
          x_c = (c_double * len(x))(*x)
          y_c = (c_double * len(y))(*y)
          z_c = (c_double * len(z))(*z)
          betaV_c = c_double(betaV)
          betaL_c = c_double(betaL)
          phase_c = c_int(phase)
          ph_c = POINTER(c_int)()
  
          self.s_sos_sound_velocity_2ph.argtypes = [POINTER( c_double ),
                                                    POINTER( c_double ),
                                                    POINTER( c_double ),
                                                    POINTER( c_double ),
                                                    POINTER( c_double ),
                                                    POINTER( c_double ),
                                                    POINTER( c_double ),
                                                    POINTER( c_int ),
                                                    POINTER( c_int )]
  
          self.s_sos_sound_velocity_2ph.restype = c_double
  
          sos = self.s_sos_sound_velocity_2ph(byref(temp_c),
                                              byref(press_c),
                                              x_c,
                                              y_c,
                                              z_c,
                                              byref(betaV_c),
                                              byref(betaL_c),
                                              byref(phase_c),
                                              ph_c)
  
          return sos
  
      #################################
      # Flash interfaces
      #################################
  
      def set_ph_tolerance(self, tol):
          """Set tolerance of isobaric-isentalpic (PH) flash
  
          Args:
              tol (float): Tolerance
          """
          tol_c = c_double(tol)
          self.s_set_ph_tolerance.argtypes = [POINTER( c_double )]
          self.s_set_ph_tolerance.restype = None
          self.s_set_ph_tolerance(byref(tol_c))
  
      def two_phase_tpflash(self,temp,press,z):
          """Do isothermal-isobaric (TP) flash
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              z (array_like): Overall molar composition
  
          Returns:
              x (ndarray): Liquid molar composition
              y (ndarray): Gas molar composition
              betaV (float): Molar gas phase fraction
              betaL (float): Molar liquid phase fraction
              phase (int): Phase identifier (iTWOPH/iLIQPH/iVAPPH)
          """
          self.activate()
          temp_c = c_double(temp)
          press_c = c_double(press)
          z_c = (c_double * len(z))(*z)
  
          x_c = (c_double * len(z))(0.0)
          y_c = (c_double * len(z))(0.0)
          betaV_c = c_double(0.0)
          betaL_c = c_double(0.0)
          phase_c = c_int(0)
  
          self.s_twophasetpflash.argtypes = [POINTER( c_double ),
                                             POINTER( c_double ),
                                             POINTER( c_double ),
                                             POINTER( c_double ),
                                             POINTER( c_double ),
                                             POINTER( c_int ),
                                             POINTER( c_double ),
                                             POINTER( c_double )]
  
          self.s_twophasetpflash.restype = None
  
          self.s_twophasetpflash(byref(temp_c),
                                 byref(press_c),
                                 z_c,
                                 byref(betaV_c),
                                 byref(betaL_c),
                                 byref(phase_c),
                                 x_c,
                                 y_c)
  
          x = np.array(x_c)
          y = np.array(y_c)
  
          return x, y, betaV_c.value, betaL_c.value, phase_c.value
  
  
      def two_phase_psflash(self,press,z,entropy,temp=None):
          """Do isentropic-isobaric (SP) flash
  
          Args:
              press (float): Pressure (Pa)
              z (array_like): Overall molar composition
              entropy (float): Specific entropy (J/mol/K)
              temp (float, optional): Initial guess for temperature (K)
  
          Returns:
              temp (float): Temperature (K)
              x (ndarray): Liquid molar composition
              y (ndarray): Gas molar composition
              betaV (float): Molar gas phase fraction
              betaL (float): Molar liquid phase fraction
              phase (int): Phase identifier (iTWOPH/iLIQPH/iVAPPH)
          """
          self.activate()
          press_c = c_double(press)
          z_c = (c_double * len(z))(*z)
          s_c = c_double(entropy)
  
          if not temp is None:
              temp_c = POINTER( c_double )(c_double(temp))
          else:
              temp_c = POINTER( c_double )(c_double(0.0))
  
          x_c = (c_double * len(z))(0.0)
          y_c = (c_double * len(z))(0.0)
          betaV_c = c_double(0.0)
          betaL_c = c_double(0.0)
          phase_c = c_int(0)
          ierr_c = c_int(0)
          self.s_psflash_twophase.argtypes = [POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_int ),
                                              POINTER( c_int )]
  
          self.s_psflash_twophase.restype = None
  
          self.s_psflash_twophase(temp_c,
                                  byref(press_c),
                                  z_c,
                                  byref(betaV_c),
                                  byref(betaL_c),
                                  x_c,
                                  y_c,
                                  byref(s_c),
                                  byref(phase_c),
                                  byref(ierr_c))
  
          if ierr_c.value != 0:
              raise Exception("PS flash calclualtion failed")
  
          x = np.array(x_c)
          y = np.array(y_c)
  
          return temp_c[0], x, y, betaV_c.value, betaL_c.value, phase_c.value
  
      def two_phase_phflash(self,press,z,enthalpy,temp=None):
          """Do isenthalpic-isobaric (HP) flash
  
          Args:
              press (float): Pressure (Pa)
              z (array_like): Overall molar composition
              enthalpy (float): Specific enthalpy (J/mol)
              temp (float, optional): Initial guess for temperature (K)
  
          Returns:
              temp (float): Temperature (K)
              x (ndarray): Liquid molar composition
              y (ndarray): Gas molar composition
              betaV (float): Molar gas phase fraction
              betaL (float): Molar liquid phase fraction
              phase (int): Phase identifier (iTWOPH/iLIQPH/iVAPPH)
          """
          self.activate()
          press_c = c_double(press)
          z_c = (c_double * len(z))(*z)
          h_c = c_double(enthalpy)
  
          if not temp is None:
              temp_c = POINTER( c_double )(c_double(temp))
          else:
              temp_c = POINTER( c_double )(c_double(0.0))
  
          x_c = (c_double * len(z))(0.0)
          y_c = (c_double * len(z))(0.0)
          betaV_c = c_double(0.0)
          betaL_c = c_double(0.0)
          phase_c = c_int(0)
          ierr_c = c_int(0)
  
          self.s_phflash_twophase.argtypes = [POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_int ),
                                              POINTER( c_int )]
  
          self.s_phflash_twophase.restype = None
  
          self.s_phflash_twophase(temp_c,
                                  byref(press_c),
                                  z_c,
                                  byref(betaV_c),
                                  byref(betaL_c),
                                  x_c,
                                  y_c,
                                  byref(h_c),
                                  byref(phase_c),
                                  byref(ierr_c))
  
          if ierr_c.value != 0:
              raise Exception("PH flash calclualtion failed")
  
          x = np.array(x_c)
          y = np.array(y_c)
  
          return temp_c[0], x, y, betaV_c.value, betaL_c.value, phase_c.value
  
      def two_phase_uvflash(self,z,specific_energy,specific_volume,temp=None,press=None):
          """Do isoenergetic-isochoric (UV) flash
  
          Args:
              press (float): Pressure (Pa)
              z (array_like): Overall molar composition
              specific_energy (float): Specific energy (J/mol)
              specific_volume (float): Specific volume (m3/mol)
              temp (float, optional): Initial guess for temperature (K)
              press (float, optional): Initial guess for pressure (Pa)
  
          Returns:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
              x (ndarray): Liquid molar composition
              y (ndarray): Gas molar composition
              betaV (float): Molar gas phase fraction
              betaL (float): Molar liquid phase fraction
              phase (int): Phase identifier (iTWOPH/iLIQPH/iVAPPH)
          """
          self.activate()
          z_c = (c_double * len(z))(*z)
          e_c = c_double(specific_energy)
          v_c = c_double(specific_volume)
  
          if not temp is None:
              temp_c = POINTER( c_double )(c_double(temp))
          else:
              temp_c = POINTER( c_double )(c_double(0.0))
  
          if not press is None:
              press_c = POINTER( c_double )(c_double(press))
          else:
              press_c = POINTER( c_double )(c_double(0.0))
  
          x_c = (c_double * len(z))(0.0)
          y_c = (c_double * len(z))(0.0)
          betaV_c = c_double(0.0)
          betaL_c = c_double(0.0)
          phase_c = c_int(0)
  
          self.s_uvflash_twophase.argtypes = [POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_int )]
  
          self.s_uvflash_twophase(temp_c,
                                  press_c,
                                  z_c,
                                  byref(betaV_c),
                                  byref(betaL_c),
                                  x_c,
                                  y_c,
                                  byref(e_c),
                                  byref(v_c),
                                  byref(phase_c))
  
          x = np.array(x_c)
          y = np.array(y_c)
  
          return temp_c[0], press_c[0], x, y, betaV_c.value, betaL_c.value, phase_c.value
  
      def guess_phase(self, temp, press, z):
          """If only one root exsist for the equation of state the phase type can be
          determined from either the psedo-critical volume or a volume ratio to the co-volume
  
          Args:
              temp (float): Temperature (K)
              press (float): Pressure (Pa)
  
          Returns:
              int: Phase int (VAPPH or LIQPH)
          """
          self.activate()
          temp_c = c_double(temp)
          press_c = c_double(press)
          z_c = (c_double * len(z))(*z)
          null_pointer = POINTER(c_double)()
          temp_comp_c = null_pointer
          press_comp_c = null_pointer
          vb_ratio_c = null_pointer
  
          self.s_guess_phase.argtypes = [POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double )]
  
          self.s_guess_phase.restype = c_int
  
          phase = self.s_guess_phase(byref(temp_c),
                                     byref(press_c),
                                     z_c,
                                     temp_comp_c,
                                     press_comp_c,
                                     vb_ratio_c)
  
          return phase
  
      #################################
      # Temperature-volume property interfaces
      #################################
  
      def pressure_tv(self, temp, volume, n, dpdt=None, dpdv=None, dpdn=None):
          """Calculate pressure given temperature, volume and mol numbers.
  
          Args:
              temp (float): Temperature (K)
              volume (float): Volume (m3)
              n (array_like): Mol numbers (mol)
              dpdt (No type, optional): Flag to activate calculation. Defaults to None.
              dpdv (No type, optional): Flag to activate calculation. Defaults to None.
              dpdn (No type, optional): Flag to activate calculation. Defaults to None.
  
          Returns:
              float: Pressure (Pa)
              Optionally pressure differentials
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(volume)
          n_c = (c_double * len(n))(*n)
  
          null_pointer = POINTER(c_double)()
          if dpdt is None:
              dpdt_c = null_pointer
          else:
              dpdt_c = POINTER(c_double)(c_double(0.0))
          if dpdv is None:
              dpdv_c = null_pointer
          else:
              dpdv_c = POINTER(c_double)(c_double(0.0))
          d2pdv2_c = null_pointer
          if dpdn is None:
              dpdn_c = null_pointer
          else:
              dpdn_c = (c_double * len(n))(0.0)
  
          recalculate_c = POINTER(c_int)(c_int(1))
  
          self.s_pressure_tv.argtypes = [POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_int )]
  
          self.s_pressure_tv.restype = c_double
  
          P = self.s_pressure_tv(byref(temp_c),
                                 byref(v_c),
                                 n_c,
                                 dpdv_c,
                                 dpdt_c,
                                 d2pdv2_c,
                                 dpdn_c,
                                 recalculate_c)
  
          return_tuple = (P, )
          if not dpdt is None:
              return_tuple += (dpdt_c[0], )
          if not dpdv is None:
              return_tuple += (dpdv_c[0], )
          if not dpdn is None:
              return_tuple += (np.array(dpdn_c), )
  
          return return_tuple
  
      def internal_energy_tv(self, temp, volume, n, dedt=None, dedv=None):
          """Calculate internal energy given temperature, volume and mol numbers.
  
          Args:
              temp (float): Temperature (K)
              volume (float): Volume (m3)
              n (array_like): Mol numbers (mol)
              dedt (No type, optional): Flag to activate calculation. Defaults to None.
              dedv (No type, optional): Flag to activate calculation. Defaults to None.
  
          Returns:
              float: Energy (J)
              Optionally energy differentials
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(volume)
          e_c = c_double(0.0)
          n_c = (c_double * len(n))(*n)
  
          null_pointer = POINTER(c_double)()
          if dedt is None:
              dedt_c = null_pointer
          else:
              dedt_c = POINTER(c_double)(c_double(0.0))
          if dedv is None:
              dedv_c = null_pointer
          else:
              dedv_c = POINTER(c_double)(c_double(0.0))
  
          recalculate_c = POINTER(c_int)(c_int(1))
  
          self.s_internal_energy_tv.argtypes = [POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_int )]
  
          self.s_internal_energy_tv.restype = None
  
          self.s_internal_energy_tv(byref(temp_c),
                                    byref(v_c),
                                    n_c,
                                    byref(e_c),
                                    dedt_c,
                                    dedv_c,
                                    recalculate_c)
  
          return_tuple = (e_c.value, )
          if not dedt is None:
              return_tuple += (dedt_c[0], )
          if not dedv is None:
              return_tuple += (dedv_c[0], )
  
          return return_tuple
  
      def entropy_tv(self, temp, volume, n, dsdt=None, dsdv=None, dsdn=None):
          """Calculate entropy given temperature, volume and mol numbers.
  
          Args:
              temp (float): Temperature (K)
              volume (float): Volume (m3)
              n (array_like): Mol numbers (mol)
              dsdt (No type, optional): Flag to activate calculation. Defaults to None.
              dsdv (No type, optional): Flag to activate calculation. Defaults to None.
              dsdn (No type, optional): Flag to activate calculation. Defaults to None.
  
          Returns:
              float: Entropy (J/K)
              Optionally entropy differentials
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(volume)
          s_c = c_double(0.0)
          n_c = (c_double * len(n))(*n)
  
          null_pointer = POINTER(c_double)()
          if dsdt is None:
              dsdt_c = null_pointer
          else:
              dsdt_c = POINTER(c_double)(c_double(0.0))
          if dsdv is None:
              dsdv_c = null_pointer
          else:
              dsdv_c = POINTER(c_double)(c_double(0.0))
          if dsdn is None:
              dsdn_c = null_pointer
          else:
              dsdn_c = (c_double * len(n))(0.0)
  
          residual_c = POINTER(c_int)(c_int(0))
  
          self.s_entropy_tv.argtypes = [POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_double ),
                                        POINTER( c_int )]
  
          self.s_entropy_tv.restype = None
  
          self.s_entropy_tv(byref(temp_c),
                            byref(v_c),
                            n_c,
                            byref(s_c),
                            dsdt_c,
                            dsdv_c,
                            dsdn_c,
                            residual_c)
  
          return_tuple = (s_c.value, )
          if not dsdt is None:
              return_tuple += (dsdt_c[0], )
          if not dsdv is None:
              return_tuple += (dsdv_c[0], )
          if not dsdn is None:
              return_tuple += (np.array(dsdn_c), )
  
          return return_tuple
  
      def enthalpy_tv(self, temp, volume, n, dhdt=None, dhdv=None, dhdn=None):
          """Calculate enthalpy given temperature, volume and mol numbers.
  
          Args:
              temp (float): Temperature (K)
              volume (float): Volume (m3)
              n (array_like): Mol numbers (mol)
              dhdt (No type, optional): Flag to activate calculation. Defaults to None.
              dhdv (No type, optional): Flag to activate calculation. Defaults to None.
              dhdn (No type, optional): Flag to activate calculation. Defaults to None.
  
          Returns:
              float: Enthalpy (J)
              Optionally enthalpy differentials
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(volume)
          h_c = c_double(0.0)
          n_c = (c_double * len(n))(*n)
  
          null_pointer = POINTER(c_double)()
          if dhdt is None:
              dhdt_c = null_pointer
          else:
              dhdt_c = POINTER(c_double)(c_double(0.0))
          if dhdv is None:
              dhdv_c = null_pointer
          else:
              dhdv_c = POINTER(c_double)(c_double(0.0))
          if dhdn is None:
              dhdn_c = null_pointer
          else:
              dhdn_c = (c_double * len(n))(0.0)
  
          residual_c = POINTER(c_int)(c_int(0))
  
          self.s_enthalpy_tv.argtypes = [POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_int )]
  
          self.s_enthalpy_tv.restype = None
  
          self.s_enthalpy_tv(byref(temp_c),
                             byref(v_c),
                             n_c,
                             byref(h_c),
                             dhdt_c,
                             dhdv_c,
                             dhdn_c,
                             residual_c)
  
          return_tuple = (h_c.value, )
          if not dhdt is None:
              return_tuple += (dhdt_c[0], )
          if not dhdv is None:
              return_tuple += (dhdv_c[0], )
          if not dhdn is None:
              return_tuple += (np.array(dhdn_c), )
  
          return return_tuple
  
      def helmholtz_tv(self, temp, volume, n, dadt=None, dadv=None):
          """Calculate Helmholtz energy given temperature, volume and mol numbers.
  
          Args:
              temp (float): Temperature (K)
              volume (float): Volume (m3)
              n (array_like): Mol numbers (mol)
              dadt (No type, optional): Flag to activate calculation. Defaults to None.
              dadv (No type, optional): Flag to activate calculation. Defaults to None.
  
          Returns:
              float: Helmholtz energy (J)
              Optionally energy differentials
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(volume)
          a_c = c_double(0.0)
          n_c = (c_double * len(n))(*n)
  
          null_pointer = POINTER(c_double)()
          if dadt is None:
              dadt_c = null_pointer
          else:
              dadt_c = POINTER(c_double)(c_double(0.0))
          if dadv is None:
              dadv_c = null_pointer
          else:
              dadv_c = POINTER(c_double)(c_double(0.0))
          d2adt2_c = null_pointer
          d2adv2_c = null_pointer
          d2advdt_c = null_pointer
          recalculate_c = POINTER(c_int)(c_int(1))
  
          self.s_helmholtz_energy.argtypes = [POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_double ),
                                              POINTER( c_int )]
  
          self.s_helmholtz_energy.restype = None
  
          self.s_helmholtz_energy(byref(temp_c),
                                  byref(v_c),
                                  n_c,
                                  byref(a_c),
                                  dadt_c,
                                  dadv_c,
                                  d2adt2_c,
                                  d2adv2_c,
                                  d2advdt_c,
                                  recalculate_c)
  
          return_tuple = (a_c.value, )
          if not dadt is None:
              return_tuple += (dadt_c[0], )
          if not dadv is None:
              return_tuple += (dadv_c[0], )
  
          return return_tuple
  
      def chemical_potential_tv(self, temp, volume, n, dmudt=None, dmudv=None, dmudn=None):
          """Calculate chemical potential given temperature, volume and mol numbers.
  
          Args:
              temp (float): Temperature (K)
              volume (float): Volume (m3)
              n (array_like): Mol numbers (mol)
              dmudt (No type, optional): Flag to activate calculation. Defaults to None.
              dmudv (No type, optional): Flag to activate calculation. Defaults to None.
              dmudn (No type, optional): Flag to activate calculation. Defaults to None.
  
          Returns:
              float: Chemical potential (J/mol)
              Optionally chemical potential differentials
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(volume)
          mu_c = (c_double * len(n))(0.0)
          n_c = (c_double * len(n))(*n)
  
          null_pointer = POINTER(c_double)()
          if dmudt is None:
              dmudt_c = null_pointer
          else:
              dmudt_c = (c_double * len(n))(0.0)
          if dmudv is None:
              dmudv_c = null_pointer
          else:
              dmudv_c = (c_double * len(n))(0.0)
          if dmudn is None:
              dmudn_c = null_pointer
          else:
              dmudn_c = (c_double * len(n)**2)(0.0)
  
          recalculate_c = POINTER(c_int)(c_int(1))
  
          self.s_chempot.argtypes = [POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_int )]
  
          self.s_chempot.restype = None
  
          self.s_chempot(byref(temp_c),
                         byref(v_c),
                         n_c,
                         mu_c,
                         dmudv_c,
                         dmudt_c,
                         dmudn_c,
                         recalculate_c)
  
          return_tuple = (np.array(mu_c), )
          if not dmudt is None:
              return_tuple += (np.array(dmudt_c), )
          if not dmudv is None:
              return_tuple += (np.array(dmudv_c), )
          if not dmudn is None:
              dmudn = np.zeros((len(n), len(n)))
              for i in range(len(n)):
                  for j in range(len(n)):
                      dmudn[i][j] = dmudn_c[i + j*len(n)]
              return_tuple += (np.array(dmudn), )
  
          return return_tuple
  
      def fugacity_tv(self, temp, volume, n, dlnphidt=None, dlnphidv=None, dlnphidn=None):
          """Calculate natural logarithm of fugacity given temperature, volume and mol numbers.
  
          Args:
              temp (float): Temperature (K)
              volume (float): Volume (m3)
              n (array_like): Mol numbers (mol)
              dlnphidt (No type, optional): Flag to activate calculation. Defaults to None.
              dlnphidv (No type, optional): Flag to activate calculation. Defaults to None.
              dlnphidn (No type, optional): Flag to activate calculation. Defaults to None.
  
          Returns:
              ndarry: Natural logarithm of fugacity
              Optionally differentials
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(volume)
          lnphi_c = (c_double * len(n))(0.0)
          n_c = (c_double * len(n))(*n)
  
          null_pointer = POINTER(c_double)()
          if dlnphidt is None:
              dlnphidt_c = null_pointer
          else:
              dlnphidt_c = (c_double * len(n))(0.0)
          if dlnphidv is None:
              dlnphidv_c = null_pointer
          else:
              dlnphidv_c = (c_double * len(n))(0.0)
          if dlnphidn is None:
              dlnphidn_c = null_pointer
          else:
              dlnphidn_c = (c_double * len(n)**2)(0.0)
  
          self.s_lnphi_tv.argtypes = [POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double )]
  
          self.s_lnphi_tv.restype = None
  
          self.s_lnphi_tv(byref(temp_c),
                          byref(v_c),
                          n_c,
                          lnphi_c,
                          dlnphidt_c,
                          dlnphidv_c,
                          dlnphidn_c)
  
          return_tuple = (np.array(lnphi_c), )
          if not dlnphidt is None:
              return_tuple += (np.array(dlnphidt_c), )
          if not dlnphidv is None:
              return_tuple += (np.array(dlnphidv_c), )
          if not dlnphidn is None:
              dlnphidn = np.zeros((len(n),len(n)))
              for i in range(len(n)):
                  for j in range(len(n)):
                      dlnphidn[i][j] = dlnphidn_c[i + j*len(n)]
              return_tuple += (dlnphidn, )
  
          return return_tuple
  
      #################################
      # Saturation interfaces
      #################################
  
      def bubble_temperature(self, press, z):
          """Calculate bubble temperature given pressure and composition
  
          Args:
              press (float): Pressure (Pa)
              z (array_like): Composition (-)
  
          Raises:
              Exception: Faild to calculate
  
          Returns:
              float: Temperature (K)
              ndarray: Incipient phase composition
          """
          self.activate()
          press_c = c_double(press)
          y_c = (c_double * len(z))(0.0)
          z_c = (c_double * len(z))(*z)
          ierr_c = c_int(0)
  
          self.s_bubble_t.argtypes = [POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_int )]
  
          self.s_bubble_t.restype = c_double
  
          temp = self.s_bubble_t(byref(press_c),
                                 z_c,
                                 y_c,
                                 byref(ierr_c))
  
          y = np.array(y_c)
          if ierr_c.value != 0:
              raise Exception("bubble_temperature calclualtion failed")
          return temp, y
  
      def bubble_pressure(self, temp, z):
          """Calculate bubble pressure given temperature and composition
  
          Args:
              temp (float): Temperature (K)
              z (array_like): Composition (-)
  
          Raises:
              Exception: Faild to calculate
  
          Returns:
              float: Pressure (Pa)
              ndarray: Incipient phase composition
          """
          self.activate()
          temp_c = c_double(temp)
          y_c = (c_double * len(z))(0.0)
          z_c = (c_double * len(z))(*z)
          ierr_c = c_int(0)
  
          self.s_bubble_p.argtypes = [POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_int )]
  
          self.s_bubble_p.restype = c_double
  
          press = self.s_bubble_p(byref(temp_c),
                                  z_c,
                                  y_c,
                                  byref(ierr_c))
  
          y = np.array(y_c)
          if ierr_c.value != 0:
              raise Exception("bubble_pressure calclualtion failed")
          return press, y
  
      def dew_temperature(self,press,z):
          """Calculate dew temperature given pressure and composition
  
          Args:
              temp (float): Pressure (Pa)
              z (float): Compositon (-)
  
          Raises:
              Exception: Not able to solve for dew point
  
          Returns:
              float : Temperature (K)
              ndarray : Incipient phase composition (-)
          """
          self.activate()
          press_c = c_double(press)
          x_c = (c_double * len(z))(0.0)
          z_c = (c_double * len(z))(*z)
          ierr_c = c_int(0)
  
          self.s_dew_t.argtypes = [POINTER( c_double ),
                                   POINTER( c_double ),
                                   POINTER( c_double ),
                                   POINTER( c_int )]
  
          self.s_dew_t.restype = c_double
  
          temp = self.s_dew_t(byref(press_c),
                              x_c,
                              z_c,
                              byref(ierr_c))
  
          x = np.array(x_c)
          if ierr_c.value != 0:
              raise Exception("dew_temperature calclualtion failed")
          return temp, x
  
      def dew_pressure(self,temp,z):
          """Calculate dew pressure given temperature and composition
  
          Args:
              temp (float): Temperature (K)
              z (float): Compositon (-)
  
          Raises:
              Exception: Not able to solve for dew point
  
          Returns:
              float : Pressure (Pa)
              ndarray : Incipient phase composition (-)
          """
          self.activate()
          temp_c = c_double(temp)
          x_c = (c_double * len(z))(0.0)
          z_c = (c_double * len(z))(*z)
          ierr_c = c_int(0)
  
          self.s_dew_p.argtypes = [POINTER( c_double ),
                                   POINTER( c_double ),
                                   POINTER( c_double ),
                                   POINTER( c_int )]
  
          self.s_dew_p.restype = c_double
  
          press = self.s_dew_p(byref(temp_c),
                               x_c,
                               z_c,
                               byref(ierr_c))
  
          x = np.array(x_c)
          if ierr_c.value != 0:
              raise Exception("bubble_pressure calclualtion failed")
          return press, x
  
      def get_envelope_twophase(self, initial_pressure, z, maximum_pressure=1.5e7,
                                minimum_temperature=None, step_size=None,
                                calc_v=False):
          """Get the phase-envelope
  
          Args:
              initial_pressure (float): Start mapping form dew point at initial pressure (Pa).
              z (array_like): Composition (-)
              maximum_pressure (float , optional): Exit on maximum pressure (Pa). Defaults to 1.5e7.
              minimum_temperature (float , optional): Exit on minimum pressure (Pa). Defaults to None.
              step_size (float , optional): Tune step size of envelope trace. Defaults to None.
              calc_v (bool, optional): Calculate specifc volume of saturated phase? Defaults to False
          Returns:
              ndarray: Temperature values (K)
              ndarray: Pressure values (Pa)
              ndarray (optional): Specific volume (m3/mol)
          """
          self.activate()
          nmax = 1000
          z_c = (c_double * len(z))(*z)
          temp_c = c_double(0.0)
          press_c = c_double(initial_pressure)
          spec_c = c_int(1)
          beta_in_c = c_double(1.0)
          max_press_c = c_double(maximum_pressure)
          nmax_c = c_int(nmax)
          Ta_c = (c_double * nmax)(0.0)
          Pa_c = (c_double * nmax)(0.0)
          Ki_c = (c_double * (nmax*len(z)))(0.0)
          beta_c = (c_double * nmax)(0.0)
          n_c = c_int(0)
          null_pointer = POINTER(c_double)()
          criconden_c = null_pointer
          crit_c = null_pointer
          if step_size is None:
              ds_c = null_pointer
          else:
              ds_c = POINTER(c_double)(c_double(step_size))
          exitOnTriplePoint_c = POINTER(c_int)()
          if minimum_temperature is None:
              tme_c = null_pointer
          else:
              tme_c = POINTER(c_double)(c_double(minimum_temperature))
  
          self.s_envelope_plot.argtypes = [POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_int ),
                                           POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_int ),
                                           POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_int ),
                                           POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_double ),
                                           POINTER( c_int ),
                                           POINTER( c_double )]
  
          self.s_envelope_plot.restype = None
  
          self.s_envelope_plot(z_c,
                               byref(temp_c),
                               byref(press_c),
                               byref(spec_c),
                               byref(beta_in_c),
                               byref(max_press_c),
                               byref(nmax_c),
                               Ta_c,
                               Pa_c,
                               Ki_c,
                               beta_c,
                               byref(n_c),
                               criconden_c,
                               crit_c,
                               ds_c,
                               exitOnTriplePoint_c,
                               tme_c)
  
          t_vals = np.array(Ta_c[0:n_c.value])
          p_vals = np.array(Pa_c[0:n_c.value])
  
          return_tuple = (t_vals, p_vals)
  
          if calc_v:
              # Special treatment for single phase
              if np.amax(z) == 1:
                  t_vals_single = np.zeros(2*n_c.value)
                  p_vals_single = np.zeros(2*n_c.value)
                  v_vals_single = np.zeros_like(t_vals_single)
                  for i in range(n_c.value):
                      t_vals_single[i] = t_vals[i]
                      t_vals_single[-i-1] = t_vals[i]
                      p_vals_single[i] = p_vals[i]
                      p_vals_single[-i-1] = p_vals[i]
                      v_vals_single[i], = self.specific_volume(t_vals[i], p_vals[i], z, self.VAPPH)
                      v_vals_single[-i-1], = self.specific_volume(t_vals[i], p_vals[i], z, self.LIQPH)
                  return_tuple = (t_vals_single, p_vals_single, v_vals_single)
              else:
                  v_vals = np.zeros_like(t_vals)
                  for i in range(n_c.value):
                      if beta_c[i] > 0.5:
                          phase = self.VAPPH
                      else:
                          phase = self.LIQPH
                      v_vals[i], = self.specific_volume(t_vals[i], p_vals[i], z, phase)
                  return_tuple += (v_vals, )
  
          return return_tuple
  
      def get_binary_pxy(self,
                         temp,
                         maximum_pressure=1.5e7,
                         minimum_pressure=1.0e5,
                         maximum_dz = 0.003,
                         maximum_dlns=0.01):
          """Calculate binary three phase envelope
  
          Args:
              temp (float): Temperature (K)
              maximum_pressure (float, optional): Exit on maximum pressure (Pa). Defaults to 1.5e7.
              minimum_pressure (float, optional): Exit on minimum pressure (Pa). Defaults to 1.0e5.
              maximum_dz (float, optional): [description]. Defaults to 0.003.
              maximum_dlns (float, optional): [description]. Defaults to 0.01.
  
          Returns:
              tuple of arrays: LLE, L1VE, L2VE
          """
          # Redefinition of module parameter:
          self.activate()
          nmax = 10000
          #c_int.in_dll(self.tp, self.get_export_name("binaryplot", "maxpoints")).value
  
          temp_c = c_double(temp)
          min_temp_c = c_double(0.0)
          ispec_c = c_int(1)
          press_c = c_double(0.0)
          max_press_c = c_double(maximum_pressure)
          min_press_c = c_double(minimum_pressure)
          dz_max_c = c_double(maximum_dz)
          dlns_max_c = c_double(maximum_dlns)
          filename = "binaryVLLE.dat"
          filename_c = c_char_p(filename.encode('ascii'))
          filename_len = c_len_type(len(filename))
          res_c = (c_double * (nmax*9))(0.0)
          nres_c = (c_int * 3)(0)
          wsf_c = c_int(1)
  
          self.s_binary_plot.argtypes = [POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_int ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_char_p ),
                                         POINTER( c_double ),
                                         POINTER( c_double ),
                                         POINTER( c_int ),
                                         POINTER( c_int ),
                                         POINTER( c_double ),
                                         c_len_type]
  
          self.s_binary_plot.restype = None
  
          self.s_binary_plot(byref(temp_c),
                             byref(press_c),
                             byref(ispec_c),
                             byref(min_temp_c),
                             byref(max_press_c),
                             byref(dz_max_c),
                             filename_c,
                             byref(dlns_max_c),
                             res_c,
                             nres_c,
                             byref(wsf_c),
                             byref(min_press_c),
                             filename_len)
  
          nLLE = nres_c[0]
          nL1VE = nres_c[1]
          nL2VE = nres_c[2]
  
          if nLLE > 0:
              xLLE = np.zeros(nLLE)
              wLLE = np.zeros(nLLE)
              pLLE = np.zeros(nLLE)
              for i in range(nLLE):
                  xLLE[i] = res_c[i*9]
                  wLLE[i] = res_c[i*9+1]
                  pLLE[i] = res_c[i*9+2]
              LLE = (xLLE, wLLE, pLLE)
          else:
              LLE = (None, None, None)
  
          if nL1VE > 0:
              xL1VE = np.zeros(nL1VE)
              wL1VE = np.zeros(nL1VE)
              pL1VE = np.zeros(nL1VE)
              for i in range(nL1VE):
                  xL1VE[i] = res_c[i*9+3]
                  wL1VE[i] = res_c[i*9+4]
                  pL1VE[i] = res_c[i*9+5]
              L1VE = (xL1VE, wL1VE, pL1VE)
          else:
              L1VE = (None, None, None)
  
          if nL2VE > 0:
              xL2VE = np.zeros(nL2VE)
              wL2VE = np.zeros(nL2VE)
              pL2VE = np.zeros(nL2VE)
              for i in range(nL2VE):
                  xL2VE[i] = res_c[i*9+6]
                  wL2VE[i] = res_c[i*9+7]
                  pL2VE[i] = res_c[i*9+8]
              L2VE = (xL2VE, wL2VE, pL2VE)
          else:
              L2VE = (None, None, None)
  
          return LLE, L1VE, L2VE
  
      def get_bp_term(self,
                      i_term):
          """Get error description for binary plot error
  
          Args:
              i_term (int): binary plot error identifyer
  
          Returns:
              str: Error message
          """
          message_len = 50
          message_c = c_char_p(b" " * message_len)
          message_len = c_len_type(message_len)
          i_term_c = c_int(i_term)
  
          self.s_get_bp_term.argtypes = [POINTER( c_int ),
                                         c_char_p,
                                         c_len_type]
  
          self.s_get_bp_term.restype = None
  
          self.s_get_bp_term(byref(i_term_c),
                             message_c,
                             message_len)
          message = message_c.value.decode('ascii')
          return message
  
      def global_binary_plot(self,
                             maximum_pressure=1.5e7,
                             minimum_pressure=1.0e5,
                             minimum_temperature=150.0,
                             maximum_temperature=500.0,
                             include_azeotropes=False):
          """Calculate global binary phase envelope
  
          Args:
              maximum_pressure (float, optional): Exit on maximum pressure (Pa). Defaults to 1.5e7.
              minimum_pressure (float, optional): Exit on minimum pressure (Pa). Defaults to 1.0e5.
              minimum_temperature (float, optional): Terminate phase line traceing at minimum temperature. Defaults to 150.0 K.
              maximum_temperature (float, optional): Terminate phase line traceing at maximum temperature. Defaults to 500.0 K.
              include_azeotropes (bool, optional): Include azeotropic lines. Defaults to False.
  
          Returns:
              tuple of arrays
          """
          self.activate()
          max_press_c = c_double(maximum_pressure)
          min_press_c = c_double(minimum_pressure)
          max_temp_c = c_double(maximum_temperature)
          min_temp_c = c_double(minimum_temperature)
          az_bool_c = c_int(1 if include_azeotropes else 0)
          filename = "global_binary.dat"
          filename_c = c_char_p(filename.encode('ascii'))
          filename_len = c_len_type(len(filename))
          i_term_c = c_int(0)
  
          self.s_global_binary_plot.argtypes = [POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                c_char_p,
                                                POINTER( c_int ),
                                                POINTER( c_double ),
                                                POINTER( c_int ),
                                                c_len_type]
  
          self.s_global_binary_plot.restype = None
  
          self.s_global_binary_plot(byref(min_press_c),
                                    byref(max_press_c),
                                    byref(min_temp_c),
                                    filename_c,
                                    byref(i_term_c),
                                    byref(max_temp_c),
                                    byref(az_bool_c),
                                    filename_len)
  
          if not i_term_c.value == 0:
              message = self.get_bp_term(i_term_c.value)
              print(message)
  
          # Load file with filename and read into arrays
          return plotutils.get_globa_binary_data(filename)
  
      def solid_envelope_plot(self, initial_pressure, z, maximum_pressure=1.5e7,
                              minimum_temperature=170.0, calc_esv=False):
          """Calculate phase envelope including solid lines
  
          Args:
              initial_pressure (float): Start mapping from initial pressure (Pa).
              z (array_like): Composition (-)
              maximum_pressure (float , optional): Exit on maximum pressure (Pa). Defaults to 1.5e7.
              calc_esv (bool, optional): Calculate specifc volume of saturated phase? Defaults to False
  
          Returns:
              tuple of arrays
          """
          self.activate()
          z_c = (c_double * len(z))(*z)
          temp_c = c_double(0.0)
          press_c = c_double(initial_pressure)
          max_press_c = c_double(maximum_pressure)
          filename = "solid_envelope.dat"
          filename_c = c_char_p(filename.encode('ascii'))
          filename_len = c_len_type(len(filename))
          i_spec_c = c_int(1)
          esv_bool_c = c_int(1 if calc_esv else 0)
  
          min_t = self.get_tmin()
          self.set_tmin(minimum_temperature)
  
          self.s_solid_envelope_plot.argtypes = [POINTER( c_double ),
                                                 POINTER( c_double ),
                                                 POINTER( c_double ),
                                                 POINTER( c_int ),
                                                 POINTER( c_double ),
                                                 c_char_p,
                                                 POINTER( c_int ),
                                                 c_len_type]
  
          self.s_solid_envelope_plot.restype = None
  
          self.s_solid_envelope_plot(z_c,
                                     byref(temp_c),
                                     byref(press_c),
                                     byref(i_spec_c),
                                     byref(max_press_c),
                                     filename_c,
                                     byref(esv_bool_c),
                                     filename_len)
  
          self.set_tmin(min_t)
  
          #if .not. i_term_c.value == 0:
          #    message = self.get_bp_term(iTerm)
          #    print(message)
  
          # Load file with filename and read into lists....
          return plotutils.get_solid_envelope_data(filename)
  
      def get_isotherm(self,
                       temp,
                       z,
                       minimum_pressure=1.0e5,
                       maximum_pressure=1.5e7,
                       nmax=100):
          """Get iso-therm at specified temperature
  
          Args:
              temp (float): Temperature (K)
              z (array_like): Composition (-)
              minimum_pressure (float, optional): Map to minimum pressure. Defaults to 1.0e5. (Pa)
              maximum_pressure (float, optional): Map to maximum pressure. Defaults to 1.5e7. (Pa)
              nmax (int, optional): Maximum number of points on iso-therm. Defaults to 100.
  
          Returns:
             Multiple numpy arrays.
          """
          self.activate()
          temp_c = c_double(temp)
          minimum_pressure_c = c_double(minimum_pressure)
          maximum_pressure_c = c_double(maximum_pressure)
          z_c = (c_double * len(z))(*z)
          va_c = (c_double * nmax)(0.0)
          pa_c = (c_double * nmax)(0.0)
          sa_c = (c_double * nmax)(0.0)
          ha_c = (c_double * nmax)(0.0)
          nmax_c = c_int(nmax)
          na_c = c_int(0)
  
          self.s_isotherm.argtypes = [POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_int ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_double ),
                                      POINTER( c_int )]
  
          self.s_isotherm.restype = None
  
          self.s_isotherm(byref(temp_c),
                          byref(minimum_pressure_c),
                          byref(maximum_pressure_c),
                          z_c,
                          byref(nmax_c),
                          pa_c,
                          va_c,
                          sa_c,
                          ha_c,
                          byref(na_c))
  
          p_vals = np.array(pa_c[0:na_c.value])
          v_vals = np.array(va_c[0:na_c.value])
          s_vals = np.array(sa_c[0:na_c.value])
          h_vals = np.array(ha_c[0:na_c.value])
  
          return p_vals, v_vals, s_vals, h_vals
  
      def get_isobar(self,
                     press,
                     z,
                     minimum_temperature=200.0,
                     maximum_temperature=500.0,
                     nmax=100):
          """Get isobar at specified pressure.
  
          Args:
              press (float): Pressure (Pa)
              z (array_like): Composition (-)
              minimum_temperature (float, optional): Minimum temperature. Defaults to 200.0. (K)
              maximum_temperature (float, optional): Maximum temperature. Defaults to 500.0. (K)
              nmax (int, optional): Maximum number of points on iso-bar. Defaults to 100.
  
          Returns:
              Multiple numpy arrays.
          """
          self.activate()
          press_c = c_double(press)
          minimum_temperature_c = c_double(minimum_temperature)
          maximum_temperature_c = c_double(maximum_temperature)
          z_c = (c_double * len(z))(*z)
          va_c = (c_double * nmax)(0.0)
          ta_c = (c_double * nmax)(0.0)
          sa_c = (c_double * nmax)(0.0)
          ha_c = (c_double * nmax)(0.0)
          nmax_c = c_int(nmax)
          na_c = c_int(0)
  
          self.s_isobar.argtypes = [POINTER( c_double ),
                                    POINTER( c_double ),
                                    POINTER( c_double ),
                                    POINTER( c_double ),
                                    POINTER( c_int ),
                                    POINTER( c_double ),
                                    POINTER( c_double ),
                                    POINTER( c_double ),
                                    POINTER( c_double ),
                                    POINTER( c_int )]
  
          self.s_isobar.restype = None
  
          self.s_isobar(byref(press_c),
                        byref(minimum_temperature_c),
                        byref(maximum_temperature_c),
                        z_c,
                        byref(nmax_c),
                        ta_c,
                        va_c,
                        sa_c,
                        ha_c,
                        byref(na_c))
  
          t_vals = np.array(ta_c[0:na_c.value])
          v_vals = np.array(va_c[0:na_c.value])
          s_vals = np.array(sa_c[0:na_c.value])
          h_vals = np.array(ha_c[0:na_c.value])
  
          return t_vals, v_vals, s_vals, h_vals
  
      def get_isenthalp(self,
                        enthalpy,
                        z,
                        minimum_pressure=1.0e5,
                        maximum_pressure=1.5e7,
                        minimum_temperature=200.0,
                        maximum_temperature=500.0,
                        nmax=100):
          """Get isenthalpy given specified enthalpy.
  
          Args:
              enthalpy (float): Enthalpy (J/mol)
              z (array_like): Composition (-)
              minimum_pressure (float, optional): Minimum pressure. Defaults to 1.0e5. (Pa)
              maximum_pressure (float, optional): Maximum pressure. Defaults to 1.5e7. (Pa)
              minimum_temperature (float, optional): Minimum temperature. Defaults to 200.0. (K)
              maximum_temperature (float, optional): Maximum temperature. Defaults to 500.0. (K)
              nmax (int, optional): Maximum number of points on isenthalp. Defaults to 100.
  
          Returns:
              Multiple numpy arrays.
          """
          self.activate()
          enthalpy_c = c_double(enthalpy)
          minimum_pressure_c = c_double(minimum_pressure)
          maximum_pressure_c = c_double(maximum_pressure)
          minimum_temperature_c = c_double(minimum_temperature)
          maximum_temperature_c = c_double(maximum_temperature)
          z_c = (c_double * len(z))(*z)
          va_c = (c_double * nmax)(0.0)
          ta_c = (c_double * nmax)(0.0)
          sa_c = (c_double * nmax)(0.0)
          pa_c = (c_double * nmax)(0.0)
          nmax_c = c_int(nmax)
          na_c = c_int(0)
  
          self.s_isenthalp.argtypes = [POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_int ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_int )]
  
          self.s_isenthalp.restype = None
  
          self.s_isenthalp(byref(enthalpy_c),
                           byref(minimum_pressure_c),
                           byref(maximum_pressure_c),
                           byref(minimum_temperature_c),
                           byref(maximum_temperature_c),
                           z_c,
                           byref(nmax_c),
                           pa_c,
                           va_c,
                           sa_c,
                           ta_c,
                           byref(na_c))
  
          t_vals = np.array(ta_c[0:na_c.value])
          v_vals = np.array(va_c[0:na_c.value])
          s_vals = np.array(sa_c[0:na_c.value])
          p_vals = np.array(pa_c[0:na_c.value])
  
          return t_vals, p_vals, v_vals, s_vals
  
      def get_isentrope(self,
                        entropy,
                        z,
                        minimum_pressure=1.0e5,
                        maximum_pressure=1.5e7,
                        minimum_temperature=200.0,
                        maximum_temperature=500.0,
                        nmax=100):
          """Get isentrope at specified entropy.
  
          Args:
              entropy (float): Entropy (J/mol/K)
              z (array_like): Composition (-)
              minimum_pressure (float, optional): Minimum pressure. Defaults to 1.0e5. (Pa)
              maximum_pressure (float, optional): Maximum pressure. Defaults to 1.5e7. (Pa)
              minimum_temperature (float, optional): Minimum temperature. Defaults to 200.0. (K)
              maximum_temperature (float, optional): Maximum temperature. Defaults to 500.0. (K)
              nmax (int, optional): Maximum number of points on isentrope. Defaults to 100.
  
          Returns:
              Multiple numpy arrays.
          """
          self.activate()
          entropy_c = c_double(entropy)
          minimum_pressure_c = c_double(minimum_pressure)
          maximum_pressure_c = c_double(maximum_pressure)
          minimum_temperature_c = c_double(minimum_temperature)
          maximum_temperature_c = c_double(maximum_temperature)
          z_c = (c_double * len(z))(*z)
          va_c = (c_double * nmax)(0.0)
          ta_c = (c_double * nmax)(0.0)
          ha_c = (c_double * nmax)(0.0)
          pa_c = (c_double * nmax)(0.0)
          nmax_c = c_int(nmax)
          na_c = c_int(0)
  
          self.s_isentrope.argtypes = [POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_int ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_double ),
                                       POINTER( c_int )]
  
          self.s_isentrope.restype = None
  
          self.s_isentrope(byref(entropy_c),
                           byref(minimum_pressure_c),
                           byref(maximum_pressure_c),
                           byref(minimum_temperature_c),
                           byref(maximum_temperature_c),
                           z_c,
                           byref(nmax_c),
                           pa_c,
                           va_c,
                           ha_c,
                           ta_c,
                           byref(na_c))
  
          t_vals = np.array(ta_c[0:na_c.value])
          v_vals = np.array(va_c[0:na_c.value])
          h_vals = np.array(ha_c[0:na_c.value])
          p_vals = np.array(pa_c[0:na_c.value])
  
          return t_vals, p_vals, v_vals, h_vals
  
      #################################
      # Stability interfaces
      #################################
  
      def critical(self, n, temp=0.0, v=0.0, tol=1.0e-7):
          """Calculate critical point in variables T and V
  
          Args:
              n (array_like): Mol numbers (mol)
              temp (float, optional): Initial guess for temperature (K). Defaults to 0.0.
              v (float, optional): Initial guess for volume (m3). Defaults to 0.0.
              tol (float, optional): Error tolerance (-). Defaults to 1.0e-8.
  
          Raises:
              Exception: Failure to solve for critcal point
  
          Returns:
              float: Temperature (K)
              float: Volume (m3)
              float: Pressure (Pa)
          """
          self.activate()
          temp_c = c_double(temp)
          v_c = c_double(v)
          n_c = (c_double * len(n))(*n)
          ierr_c = c_int(0)
          P_c = c_double(0.0)
          tol_c = c_double(tol)
          self.s_crit_tv.argtypes = [POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_double ),
                                     POINTER( c_int ),
                                     POINTER( c_double ),
                                     POINTER( c_double )]
  
          self.s_crit_tv.restype = None
  
          self.s_crit_tv(byref(temp_c),
                         byref(v_c),
                         n_c,
                         byref(ierr_c),
                         byref(tol_c),
                         byref(P_c))
  
          if ierr_c.value != 0:
              raise Exception("critical calclualtion failed")
  
          return temp_c.value, v_c.value, P_c.value
  
      #################################
      # Virial interfaces
      #################################
  
      def virial_coeffcients(self, temp, n):
          """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.
  
          Args:
              temp (float): Temperature
              n (array_like): Mol numbers (mol)
  
          Returns:
              float: B (m3/mol)
              float: C (m6/mol2)
          """
          self.activate()
          temp_c = POINTER( c_double )(c_double(temp))
          n_c = (c_double * len(n))(*n)
          B_c = c_double(0.0)
          C_c = c_double(0.0)
          self.s_virial_coeffcients.argtypes = [POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double ),
                                                POINTER( c_double )]
  
          self.s_virial_coeffcients.restype = None
  
          self.s_virial_coeffcients(temp_c,
                                    n_c,
                                    byref(B_c),
                                    byref(C_c))
  
          return B_c.value, C_c.value
  
      def second_virial_matrix(self, temp):
          """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.
  
          Args:
              temp (float): Temperature (K)
  
          Returns:
              ndarray: B - Second virial coefficient matrix (m3/mol)
          """
          self.activate()
          temp_c = POINTER( c_double )(c_double(temp))
          bmat_c = (c_double * self.nc**2)(0.0)
  
          self.s_second_virial_matrix.argtypes = [POINTER( c_double ),
                                                  POINTER( c_double )]
  
          self.s_second_virial_matrix.restype = None
  
          self.s_second_virial_matrix(temp_c, bmat_c)
  
          bmat = np.zeros((self.nc,self.nc))
          for i in range(self.nc):
              for j in range(self.nc):
                  bmat[i][j] = bmat_c[i+j*self.nc]
  
          return bmat
  
      def binary_third_virial_matrix(self, temp):
          """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
  
          Args:
              temp (float): Temperature (K)
  
          Returns:
              ndarray: C - Third virial coefficient matrix (m6/mol2)
          """
          self.activate()
          assert self.nc == 2
          temp_c = POINTER( c_double )(c_double(temp))
          cmat_c = (c_double * self.nc**2)(0.0)
  
          self.s_binary_third_virial_matrix.argtypes = [POINTER( c_double ),
                                                        POINTER( c_double )]
  
          self.s_binary_third_virial_matrix.restype = None
  
          self.s_binary_third_virial_matrix(temp_c, cmat_c)
  
          cmat = np.zeros((self.nc,self.nc))
          for i in range(self.nc):
              for j in range(self.nc):
                  cmat[i][j] = cmat_c[i+j*self.nc]
  
          return cmat
  
apis: numpy.array numpy.amax numpy.zeros_like numpy.zeros