Repository: thermotools/thermopack (language: python)

Issue:
Inital prototype / proof-of-concept for PyBind11/Fortwrap wrapping
"Reopening" PR #61 here so that we can have some nice repo structure :)


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index bac9d05f..ae4e19d5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,6 +2,7 @@ syntax: glob
 
 objects
 bin
+build/
 run*
 doc/dox.out
 doc/doxygen/html
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 00000000..356d66b8
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,6 @@
+[submodule "external/pFUnit"]
+	path = external/pFUnit
+	url = https://github.com/thermotools/pFUnit.git
+[submodule "external/lapack"]
+	path = external/lapack
+	url = https://github.com/Reference-LAPACK/lapack.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 00000000..86784bf3
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,74 @@
+string(ASCII 27 Esc)
+set(ColorRed "${Esc}[31m")
+set(ColorBlue "${Esc}[34m")
+
+cmake_minimum_required(VERSION 3.5)
+project(ThermoPack LANGUAGES Fortran)
+
+if(NOT MSVC)
+    execute_process(COMMAND bash -c "arch" OUTPUT_VARIABLE PROC)
+    set(tp_flags_common "-cpp -fPIC -fdefault-real-8 -fdefault-double-8 -frecursive -fopenmp -Wno-unused-function -Wno-unused-variable")
+
+    include(CheckCompilerFlag)
+    check_compiler_flag(Fortran "-arch arm64" arm64_supported)
+    if(arm64_supported)
+        set(gf_proc "-arch arm64")
+        set(gf_march "-arch arm64 -fno-expensive-optimizations")
+    else()
+        set(gf_proc "-mieee-fp")
+        set(gf_march "-march=x86-64 -msse2")
+    endif()
+
+    set(tp_debug_flags "${gf_proc} ${tp_flags_common} -g -fbounds-check -fbacktrace -ffpe-trap=invalid,zero,overflow -Wno-unused-dummy-argument -Wall")
+    set(tp_profile_flags "${tp_flags_common} -g -pg")
+    set(tp_optim_flags "${tp_flags_common} -O3 ${gf_march} -funroll-loops")
+
+    if(APPLE)
+        set(CMAKE_FORTRAN_COMPILER /opt/homebrew/bin/gfortran)
+        set(CMAKE_OSX_ARCHITECTURES "arm64" CACHE STRING "The OSX architecture")
+    endif()
+else()
+    # Build using: cmake --build . --config Release
+    set(tp_flags_common "/nologo /fpp /fpe:0 /names:lowercase /assume:underscore /fp:precise /extend-source:132 /iface:cref")
+    set(tp_optim_flags "${tp_flags_common}")
+    set(tp_debug_flags "${tp_flags_common} /check:bounds /traceback")
+endif()
+
+# Check if the build type is not set and set it to Release
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE)
+endif()
+
+# Set FORTRAN compile flags for thermopack and lapack
+set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} ${tp_optim_flags}")
+set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} ${tp_debug_flags}")
+
+if(MSVC)
+    if(POLICY CMP0077)
+        set(CMAKE_POLICY_DEFAULT_CMP0077 NEW)
+    endif()
+    # Disable single/complex/complex16 as it is not used by thermopack
+    set(BUILD_SINGLE OFF)
+    set(BUILD_COMPLEX OFF)
+    set(BUILD_COMPLEX16 OFF)
+    add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/lapack)
+    # Add flags not comatible with lapack
+    set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} /real-size:64 ")
+    set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} /real-size:64 /warn:all,noexternal /check:all,noarg_temp_created,nopointers")
+endif(MSVC)
+
+add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/src "thermopack")
+
+option(test "Build test suite" OFF)
+if(test)
+    cmake_policy(SET CMP0074 NEW)
+    find_package(PFUNIT QUIET)
+    if (NOT PFUNIT_FOUND)
+        message("${ColorRed}Did not find pFUnit. If you have installed pFUnit, make sure that PFUNIT_DIR is set. Falling back to building pFUnit from source...")
+        add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/pFUnit)
+    else()
+        message(STATUS "Found pFUnit: ${PFUNIT_DIR}")
+    endif()
+    message(STATUS "Creating run_unittests target")
+    add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/unittests)
+endif()
diff --git a/addon/cppExamples/.gitignore b/addon/cppExamples/.gitignore
new file mode 100644
index 00000000..567609b1
--- /dev/null
+++ b/addon/cppExamples/.gitignore
@@ -0,0 +1,1 @@
+build/
diff --git a/addon/cppExamples/CMakeLists.txt b/addon/cppExamples/CMakeLists.txt
new file mode 100644
index 00000000..f5bb2d02
--- /dev/null
+++ b/addon/cppExamples/CMakeLists.txt
@@ -0,0 +1,53 @@
+# Example for Compiling and linking C++ executables using thermopack.
+# This example assumes that libthermopack.dylib has already been built and installed with `make install`
+# The environment variable THERMOPACK_DIR should point to the root directory of ThermoPack (where thermopack-config.cmake is found)
+# this can be set with `export THERMOPACK_DIR=<path/to/thermopack>`.
+#
+# Note: On macOS (arm64), you might need to
+#   export CC=/opt/homebrew/bin/gcc-13
+#   export CXX=/opt/homebrew/bin/g++-13
+# before running cmake, in order to force compilation for arm64.
+
+project(cppExamples)
+cmake_minimum_required(VERSION 3.2)
+
+string(ASCII 27 Esc)
+set(ColourRed "${Esc}[31m")
+set(ColourBlue "${Esc}[34m")
+
+set(THERMOPACK_DIR "${CMAKE_CURRENT_SOURCE_DIR}/../..")
+
+# Searches for thermopack-config.cmake in the directory ${THERMOPACK_DIR}, which can be set with `export THERMOPACK_DIR=<path/to/thermopack>`
+# Once thermopack is found, it provides some convenience variables shown below, as well as the exported target `thermopack` which can be used for linking directly.
+find_package(THERMOPACK REQUIRED)
+if (NOT THERMOPACK_INSTALLED)
+    message("${ColourRed}Thermopack was found but is not installed. Exiting ...")
+    return()
+endif()
+
+message("${ColourBlue}Thermopack was found at : ${THERMOPACK_ROOT}")
+message("${ColourBlue}Using libthermopack at : ${THERMOPACK_LIB}")
+message("${ColourBlue}Using cppThermopack headers at : ${THERMOPACK_INCLUDE}")
+message("${ColourBlue}Linking thermopack using exported target thermopack")
+
+set(CMAKE_BUILD_TYPE "Release")
+set(CMAKE_CXX_FLAGS "-Wall -Wextra -Wfatal-errors -std=c++17")
+
+set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+set(CMAKE_CXX_FLAGS_DEBUG "-g -O0")
+
+
+add_executable(basic basic.cpp)
+target_link_libraries(basic PUBLIC thermopack)
+
+add_executable(tp_properties tp_properties.cpp)
+target_link_libraries(tp_properties PUBLIC thermopack)
+
+add_executable(tv_properties tv_properties.cpp)
+target_link_libraries(tv_properties PUBLIC thermopack)
+
+add_executable(flashes flashes.cpp)
+target_link_libraries(flashes PUBLIC thermopack)
+
+add_executable(saturation saturation.cpp)
+target_link_libraries(saturation PUBLIC thermopack)
diff --git a/addon/cppExamples/basic.cpp b/addon/cppExamples/basic.cpp
new file mode 100644
index 00000000..d5512f63
--- /dev/null
+++ b/addon/cppExamples/basic.cpp
@@ -0,0 +1,44 @@
+#include <iostream>
+#include <vector>
+#include <cppThermopack/cubic.h>
+#include <cppThermopack/saftvrmie.h>
+
+// Note: The structs `Property` and `VectorProperty` are found in utils.h, which is implicitly included in all headers
+// through thermo.h
+// Also note: Order of bools for computing differentials follows the same convention as pycThermoPack (i.e. order of bools is equal to order of input-parameters)
+// In general: The method signatures are kept consistent with pycThermopack, such that the pycThermopack docs can serve as a coarse lookup for the methods in the cpp wrapper
+//             methods that return a single value (i.e. pressure_tv, specific_volume, etc.) will return a `Property`, which has an implicit conversion to `double`, while
+//             methods that return several values (i.e. chemical_potential_tv) will return a `VectorProperty`, which has an implicit conversion to std::vector<double>
+
+int main(){
+    PengRobinson pr("N2,O2"); // Initialise Peng-Robinson EoS
+    double T = 300.;
+    double V = 1e3;
+    vector1d n = {1., 2.};
+
+    double pres = pr.pressure_tv(T, V, n); // Compute pressure (implicit conversion to double)
+    Property p = pr.pressure_tv(T, V, n, false, true, true); // Computing property and differentials
+    std::cout << "P : \n" << p << std::endl;
+    std::cout << "Pressures are equal : " << ((pres == p.value()) ? "true" : "false") << std::endl;
+    double dpdv = p.dv(); // Accessing differentials
+    std::cout << "Accessed dpdv : " << dpdv << std::endl;
+
+    std::vector<double> chem_pot = pr.chemical_potential_tv(T, V, n); // Compute chemical potential (implicit conversion to std::vector)
+    VectorProperty mu = pr.chemical_potential_tv(T, V, n, true, false, true); // Computing vector-valued property and differentials
+    std::cout << "mu :\n" << mu;
+    std::vector<double> dmudt = mu.dt(); // Accessing the differentials
+    std::vector<std::vector<double>> dmudn = mu.dn();
+
+    Saftvrmie svrm("AR,KR"); // Initialise SAFT-VR Mie EoS
+    V = 1e-3;
+    std::cout << "SVRM chemical potential : \n" << svrm.chemical_potential_tv(T, V, n, false, true) << std::endl;
+
+
+    Saftvrmie svrm_c1("C1");
+    p = 1e6;
+    T = 120.;
+    double vl = svrm_c1.specific_volume(T, p, {1}, Phase::liq);
+    double vg = svrm_c1.specific_volume(T, p, {1}, Phase::vap);
+    std::cout << "At T = " << T << " K, p = " << p / 1e5 << " bar, specific volume of methane is : " << vl << " / " << vg << std::endl;
+    return 0;
+}
\ No newline at end of file
diff --git a/addon/cppExamples/flashes.cpp b/addon/cppExamples/flashes.cpp
new file mode 100644
index 00000000..ca6ae66c
--- /dev/null
+++ b/addon/cppExamples/flashes.cpp
@@ -0,0 +1,40 @@
+#include <cppThermopack/cubic.h>
+#include <vector>
+#include <iostream>
+
+int main(){
+    SoaveRedlichKwong srk("PROP1OL,NC6");
+    std::cout << "Demonstrating flashes with SRK ...\n";
+    
+    double T1{348.}, p1{1e5};
+    std::vector<double> z = {0.4, 0.6};
+    FlashResult flsh = srk.two_phase_tpflash(T1, p1, z);
+    std::cout << "Initial state is \n" << flsh;
+
+    double h = srk.enthalpy(T1, p1, flsh.x, srk.LIQPH) * flsh.betaL 
+                + srk.enthalpy(T1, p1, flsh.y, srk.VAPPH) * flsh.betaV;
+    std::cout << "Molar enthalpy at (T, p) = (" << T1 << ", " << p1 << ") : " << h << " J / mol\n";
+    
+    double p2 = 5e3;
+    std::cout << "\nIsenthalpic decompression to " << p2 / 1e5 << " bar ...\n";
+    flsh = srk.two_phase_phflash(p2, h, z);
+    std::cout << "After isenthalpic decompression: \n" << flsh;
+
+    double s = srk.entropy(flsh.T, flsh.p, flsh.x, srk.LIQPH) * flsh.betaL
+                + srk.entropy(flsh.T, flsh.p, flsh.y, srk.VAPPH) * flsh.betaV;
+    std::cout << "Molar entropy after decompression : " << s << " J / mol K\n";
+    std::cout << "\nIsentropic compression to 2 bar ...\n";
+    flsh = srk.two_phase_psflash(2e5, s, z);
+    std::cout << "After isentropic compression : \n" << flsh;
+    
+    double v = srk.specific_volume(flsh.T, flsh.p, flsh.x, srk.LIQPH) * flsh.betaL
+                + srk.specific_volume(flsh.T, flsh.p, flsh.y, srk.VAPPH) * flsh.betaV;
+    double u = srk.internal_energy_tv(flsh.T, v, flsh.z);
+    std::cout << "Internal energy after isentropic compression : " << u << " J / mol\n";
+    
+    double q = 30.; // J / mol, Heat to add
+    std::cout << "\nIsochoric heating with " << q << " J / mol heat ...\n";
+    flsh = srk.two_phase_uvflash(u + q, v, z);
+    std::cout << "After isochoric heating : \n" << flsh;
+    return 0;
+}
\ No newline at end of file
diff --git a/addon/cppExamples/saturation.cpp b/addon/cppExamples/saturation.cpp
new file mode 100644
index 00000000..e54dcf7d
--- /dev/null
+++ b/addon/cppExamples/saturation.cpp
@@ -0,0 +1,27 @@
+#include "cppThermopack/thermo.h"
+#include "cppThermopack/pcsaft.h"
+#include <vector>
+#include <iostream>
+
+int main(){
+    Pcsaft pcs("C3,NC6,ETOH"); // Initialize pcsaft for propane/hexane/ethanol mixture
+
+    std::vector<double> z = {0.2, 0.5, 0.3}; // Total composition
+
+    std::vector<double> crit_tvp = pcs.critical(z);
+
+    double T{0.6 * crit_tvp[0]}; // Ensuring that we compute at some sub-critical temperature (0.6 * Tc)
+    std::cout << "Computing dew- and bubble pressure at " << T << " K:\n";
+    double p_dew = pcs.dew_pressure(T, z).p;
+    double p_bub = pcs.bubble_pressure(T, z).p;
+
+    std::cout << pcs.dew_pressure(T, z) << "\n";
+    std::cout << pcs.bubble_pressure(T, z) << "\n";
+
+    double p = 0.5 * (p_dew + p_bub);
+    std::cout << "\nComputing dew- and bubble temperature at " << p / 1e5 << " bar\n";
+    std::cout << pcs.dew_temperature(p, z) << "\n";
+    std::cout << pcs.bubble_temperature(p, z) << "\n";
+
+    return 0;
+}
\ No newline at end of file
diff --git a/addon/cppExamples/tp_properties.cpp b/addon/cppExamples/tp_properties.cpp
new file mode 100644
index 00000000..27f445f6
--- /dev/null
+++ b/addon/cppExamples/tp_properties.cpp
@@ -0,0 +1,23 @@
+#include <cppThermopack/saftvrmie.h>
+#include <iostream>
+#include <vector>
+
+int main(){
+    std::string comps = "C3,C1";
+    Saftvrmie svrm(comps);
+    std::cout << "Initialized SAFT-VR Mie for " << comps << "\n";
+    double T = 300;
+    double p = 1e6;
+    std::vector<double> z = {0.7, 0.3};
+
+    double h = svrm.enthalpy(T, p, z, svrm.VAPPH);
+    std::cout << "Vapour phase specific enthalpy : " << h << " J / mol\n";
+    std::cout << "Liquid phase specific enthalpy : " << static_cast<double>(svrm.enthalpy(T, p, z, svrm.LIQPH)) << " J / mol\n";
+
+    Property s_vap = svrm.entropy(T, p, z, svrm.VAPPH, false, false, true);
+    std::cout << "Vapour phase partial molar entropies : " << s_vap.dn()[0] << ", " << s_vap.dn()[0] << " J / mol^2 K\n";
+    Property s_liq = svrm.entropy(T, p, z, svrm.LIQPH, false, false, true);
+    std::cout << "Liquid phase partial molar entropies : " << s_liq.dn()[0] << ", " << s_liq.dn()[1] << " J / mol^2 K\n";
+
+    return 0;
+}
\ No newline at end of file
diff --git a/addon/cppExamples/tv_properties.cpp b/addon/cppExamples/tv_properties.cpp
new file mode 100644
index 00000000..adaa2067
--- /dev/null
+++ b/addon/cppExamples/tv_properties.cpp
@@ -0,0 +1,59 @@
+#include <cppThermopack/pcsaft.h>
+#include <iostream>
+#include <vector>
+
+// A couple convenience methods
+inline std::vector<double> operator*(double d, const std::vector<double>& v){
+    std::vector<double> v_out{v};
+    for (double& vi : v_out) vi *= d;
+    return v_out;
+}
+
+inline std::vector<double> operator*(const std::vector<double>& v, double d){
+    return d * v;
+}
+
+inline std::ostream& operator<<(std::ostream& os, const std::vector<double> v){
+    for (auto d : v) os << d << ", ";
+    return os;
+}
+
+// This is the example
+int main(){
+    std::string comps = "C1,ETOH,PROP1OL";
+    Pcsaft pcs(comps);
+    std::cout << "Initialized PC-SAFT for components " << comps << "\n";
+
+    double T{400.}, V{0.3};
+    double n_tot = 15;
+    std::vector<double> z = {0.2, 0.5, 0.5};
+    std::vector<double> n = n_tot * z;
+
+    std::cout << "Computing TV properties at T = " << T << " K, V = " << V << " m^3, n_tot = " << n_tot << " mol ... \n\n";
+
+    std::cout << "Pressure : " << pcs.pressure_tv(T, V, n) / 1e5 << " bar\n";
+    std::cout << "Internal energy : " << pcs.internal_energy_tv(T, V, n) / 1e3 << " kJ\n";
+    std::cout << "Entropy : " << pcs.entropy_tv(T, V, n) << " J / K\n";
+    std::cout << "Enthalpy : " << pcs.enthalpy_tv(T, V, n) / 1e3 << " kJ\n";
+    std::cout << "Helmholtz energy : " << pcs.helmholtz_tv(T, V, n) / 1e3 << " kJ\n";
+    std::cout << "Chemical potential : " << pcs.chemical_potential_tv(T, V, n) << " J / mol\n";
+    std::cout << "Logarithm of fugacity coeff. : " << pcs.fugacity_tv(T, V, n) << "\n";
+
+    std::cout << "\nComputing derivatives ...\n";
+    // Bools mark what derivatives to compute. 
+    // Order is dt, dv, dn
+    Property H = pcs.enthalpy_tv(T, V, n, true, true, true);
+    std::cout << "dHdT : " << H.dt() << " J / K, dHdV : " << H.dv() << " J / m^3, dHdn_i : " << H.dn() << "\n";
+    
+    Property A = pcs.helmholtz_tv(T, V, n, false, false, true);
+    std::cout << "dAdn : " << A.dn() << "\n";
+    VectorProperty mu = pcs.chemical_potential_tv(T, V, n, true, false, true);
+    std::cout << "mu   : " << mu.value() << "\n";
+    std::cout << "Chemical potential with derivatives : \n" << mu << "\n";
+
+    VectorProperty lnphi = pcs.fugacity_tv(T, V, n, false, true, true);
+    std::cout << "Logarithm of fugacity coefficients : \n" << lnphi;
+
+    std::vector<std::vector<double>> dmudn = mu.dn(); // Extracting derivatives
+    return 0;
+}
diff --git a/addon/cppThermopack/cppThermopack/cubic.h b/addon/cppThermopack/cppThermopack/cubic.h
new file mode 100644
index 00000000..1e41e07e
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/cubic.h
@@ -0,0 +1,43 @@
+#pragma once
+#include <string>
+#include "dllmacros.h"
+#include "thermo.h"
+
+extern "C" {
+    void get_export_name(eoslibinit, init_cubic)(char* comps, char* eos, char* mixing, char* alpha, char* ref,
+						int* volume_shift, size_t comp_strlen, size_t eos_strlen, size_t mixing_strlen,
+						size_t alpha_strlen, size_t ref_strlen);
+
+}
+
+class Cubic : public Thermo{
+    public:
+    Cubic(std::string comps, std::string eos,
+	    std::string mixing="vdW", std::string alpha="Classic",
+	    std::string ref="Default", bool volume_shift=false)
+	: Thermo(comps)
+    {
+	activate();
+	int vol_shift_int = volume_shift ? true_int : 0;
+	get_export_name(eoslibinit, init_cubic)(comps.data(), eos.data(), mixing.data(), alpha.data(), ref.data(),
+	    &vol_shift_int, comps.size(), eos.size(), mixing.size(), alpha.size(), ref.size());
+
+    }
+
+};
+
+class PengRobinson : public Cubic{
+    public:
+    PengRobinson(std::string comps, std::string mixing="vdW", std::string alpha="Classic",
+	    std::string ref="Default", bool volume_shift=false)
+	: Cubic(comps, "PR", mixing, alpha, ref, volume_shift)
+	{}
+};
+
+class SoaveRedlichKwong : public Cubic{
+    public:
+    SoaveRedlichKwong(std::string comps, std::string mixing="vdW", std::string alpha="Classic",
+	    std::string ref="Default", bool volume_shift=false)
+	: Cubic(comps, "SRK", mixing, alpha, ref, volume_shift)
+	{}
+};
diff --git a/addon/cppThermopack/cppThermopack/dllmacros.h b/addon/cppThermopack/cppThermopack/dllmacros.h
new file mode 100644
index 00000000..5228e7a9
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/dllmacros.h
@@ -0,0 +1,19 @@
+#pragma once
+
+typedef void* ADDRESS;
+
+#ifdef __INTEL_THERMOPACK__ // INTEL FORTRAN (x64)
+    #define MOD_PREFIX
+    #define MOD_INTERTXT _mp_
+    #define MOD_POSTFIX _
+    #define FORTRAN_POSTFIX _
+#else // GNU FORTRAN mapping
+    #define MOD_PREFIX __
+    #define MOD_INTERTXT _MOD_
+    #define MOD_POSTFIX
+    #define FORTRAN_POSTFIX _
+#endif
+
+#define _CONCAT(x, y) x##y
+#define _NEST_CONCAT(x, y) _CONCAT(x, y)
+#define get_export_name(MOD, METHOD) _NEST_CONCAT(_NEST_CONCAT(MOD_PREFIX, MOD),_NEST_CONCAT(_NEST_CONCAT(MOD_INTERTXT, METHOD), MOD_POSTFIX))
\ No newline at end of file
diff --git a/addon/cppThermopack/cppThermopack/extern.h b/addon/cppThermopack/cppThermopack/extern.h
new file mode 100644
index 00000000..8c616c4c
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/extern.h
@@ -0,0 +1,36 @@
+extern "C" {
+    int get_export_name(thermopack_var, add_eos)();
+    void get_export_name(thermopack_var, activate_model)(const int* model_idx);
+    void get_export_name(thermopack_var, delete_eos)(const int* model_idx);
+
+    void get_export_name(eos, specificvolume)(double* T, double* p, double* z, int* phase, double* v, double* dvdt, double* dvdp, double* dvdn);
+    void get_export_name(eos, zfac)(double* T, double* p, double* z, int* phase, double* zfac, double* dzdt, double* dzdp, double* dzdn);
+    void get_export_name(eos, thermo)(double* T, double* p, double* z, int* phase, double* lnfug, double* dlnfugdt, double* dlnfugdp, double* dlnfugdn, int* ophase, int* metaextremum, double* v);
+    void get_export_name(eos, enthalpy)(double* T, double* p, double* z, int* phase, double* h, double* dhdt, double* dhdp, double* dhdn, int* property_flag);
+    void get_export_name(eos, entropy)(double* T, double* p, double* z, int* phase, double* s, double* dsdt, double* dsdp, double* dsdn, int* property_flag);
+    void get_export_name(eos, ideal_enthalpy_single)(double* T, int* comp_idx, double* h_id, double* dhdt);
+    void get_export_name(eos, ideal_entropy_single)(double* T, double* p, int* comp_idx, double* s_id, double* dsdt, double* dsdp);
+    // get/set standard enthalpy/entropy/ideal cp
+    // speed_of_sound
+
+    double get_export_name(eostv, pressure)(double* T, double* V, double* n, double* dpdv, double* dpdt, double* d2pdv2, double* dpdn, int* property_flag);
+    void get_export_name(eostv, internal_energy_tv)(double* T, double* V, double* n, double* U, double* dudt, double* dudv, double* dudn, int* property_flag);
+    void get_export_name(eostv, entropy_tv)(double* T, double* V, double* n, double* S, double* dsdt, double* dsdv, double* dsdn, int* property_flag);
+    void get_export_name(eostv, enthalpy_tv)(double* T, double* V, double* n, double* S, double* dsdt, double* dsdv, double* dsdn, int* property_flag);
+    void get_export_name(eostv, free_energy_tv)(double* T, double* V, double* n, double* S, double* dsdt, double* dsdv, double* dsdn, int* property_flag);
+    void get_export_name(eostv, chemical_potential_tv)(double* T, double* V, double* n, double* mu, double* dmudt, double* dmudv, double* dmudn, int* property_flag);
+    void get_export_name(eostv, thermo_tv)(double* T, double* V, double* n, double* lnphi, double* dt, double* dv, double* dn);
+
+    void get_export_name(tp_solver, twophasetpflash)(double* T, double* p, double* z, double* betaV, double* betaL, int* phase, double* x, double* y);
+    void get_export_name(ps_solver, twophasepsflash)(double* T, double* p, double* z, double* betaV, double* betaL, double* x, double* y, double* s, int* phase, int* ierr);
+    void get_export_name(uv_solver, twophaseuvflash)(double* T, double* p, double* z, double* betaV, double* betaL, double* x, double* y, double* u, double* v, int* phase);
+    void get_export_name(ph_solver, twophasephflash)(double* T, double* p, double* z, double* betaV, double* betaL, double* x, double* y, double* h, int* phase, int* ierr);
+    // guess_phase
+
+    double get_export_name(saturation, safe_bubt)(double* p, double* z, double* y, int* ierr);
+    double get_export_name(saturation, safe_bubp)(double* T, double* z, double* y, int* ierr);
+    double get_export_name(saturation, safe_dewt)(double* p, double* z, double* x, int* ierr);
+    double get_export_name(saturation, safe_dewp)(double* T, double* z, double* x, int* ierr);
+
+    void get_export_name(critical, calccriticaltv)(double* T, double* V, double* n, int* ierr, double* tol, double* v_min, double* p);
+}
\ No newline at end of file
diff --git a/addon/cppThermopack/cppThermopack/pcsaft.h b/addon/cppThermopack/cppThermopack/pcsaft.h
new file mode 100644
index 00000000..2bcfd8ee
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/pcsaft.h
@@ -0,0 +1,31 @@
+#include "saft.h"
+
+extern "C" {
+    void get_export_name(eoslibinit, init_pcsaft)(char* comps, char* param_ref, int* simplified, int* polar, size_t comp_strlen, size_t ref_strlen);
+}
+
+class Pcsaft : public Saft{
+    public:
+    Pcsaft(std::string comps, std::string parameter_reference="Default", bool simplified=false, bool polar=false)
+        : Saft(comps)
+    {
+        activate();
+        int i_simplified = (simplified) ? true_int : 0;
+        int i_polar = (polar) ? true_int : 0;
+        get_export_name(eoslibinit, init_pcsaft)(comps.data(), parameter_reference.data(), &i_simplified, &i_polar, comps.size(), parameter_reference.size());
+    }
+};
+
+class Spcsaft : public Pcsaft{
+    public:
+    Spcsaft(std::string comps, std::string parameter_reference="Default")
+        : Pcsaft(comps, parameter_reference, true)
+    {}
+};
+
+class Pcpsaft : public Pcsaft{
+    public:
+    Pcpsaft(std::string comps, std::string parameter_reference="Default")
+        : Pcsaft(comps, parameter_reference, false, true)
+    {}
+};
\ No newline at end of file
diff --git a/addon/cppThermopack/cppThermopack/saft.h b/addon/cppThermopack/cppThermopack/saft.h
new file mode 100644
index 00000000..5ef8719b
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/saft.h
@@ -0,0 +1,8 @@
+#pragma once
+#include "thermo.h"
+
+class Saft : public Thermo {
+    public:
+    Saft(std::string& comps) : Thermo(comps){}
+
+};
\ No newline at end of file
diff --git a/addon/cppThermopack/cppThermopack/saftvrmie.h b/addon/cppThermopack/cppThermopack/saftvrmie.h
new file mode 100644
index 00000000..fc3db8f7
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/saftvrmie.h
@@ -0,0 +1,18 @@
+#pragma once
+#include <string>
+#include "dllmacros.h"
+#include "saft.h"
+
+extern "C" {
+    void get_export_name(eoslibinit, init_saftvrmie)(char* comps, char* parameter_ref, size_t comps_strlen, size_t ref_strlen);
+}
+
+class Saftvrmie : public Saft {
+    public:
+    Saftvrmie(std::string comps, std::string parameter_reference="Default")
+	: Saft(comps)
+    {
+	activate();
+	get_export_name(eoslibinit, init_saftvrmie)(comps.data(), parameter_reference.data(), comps.size(), parameter_reference.size());
+    }
+};
diff --git a/addon/cppThermopack/cppThermopack/thermo.h b/addon/cppThermopack/cppThermopack/thermo.h
new file mode 100644
index 00000000..a286987c
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/thermo.h
@@ -0,0 +1,238 @@
+#pragma once
+#include "dllmacros.h"
+#include "extern.h"
+#include "utils.h"
+#include <utility>
+#include <vector>
+#include <algorithm>
+#include <iostream>
+
+using vector1d = std::vector<double>;
+using vector2d = std::vector<vector1d>;
+
+class Thermo{
+    public:
+
+    Thermo(std::string& comps) : true_int{1}, nc{static_cast<size_t>(std::count_if(comps.begin(), comps.end(), [](char c) {return c == ',';})) + 1},
+	    model_idx{get_export_name(thermopack_var, add_eos)()}
+	{}
+
+    Thermo(Thermo&& other) : true_int{other.true_int}, nc{other.nc}, model_idx{std::exchange(other.model_idx, 0)} {}
+    Thermo(const Thermo& other) = delete; // Suggested implementation if needed : true_int{other.true_int}, nc{other.nc}, model_idx{get_export_name(thermopack_var, add_eos)()}{}
+
+    ~Thermo(){
+	// To handle move semantics (which we need), the move constructor sets the moved-from objects model_idx to
+	// zero, so that we can check in the destructor whether the object has been moved from before deleting
+	// the underlying memory with delete_eos.
+	if (model_idx > 0) get_export_name(thermopack_var, delete_eos)(&model_idx);
+    }
+
+    int TWOPH = Phase::two;
+    int LIQPH = Phase::liq;
+    int VAPPH = Phase::vap;
+
+/************************************************************************/
+// ------------------------ TP-Property Interfaces -------------------- //
+/************************************************************************/
+
+    Property specific_volume(double T, double p, vector1d z, int phase, bool dvdt=false, bool dvdp=false, bool dvdn=false) const {
+	activate();
+	Property v(nc, dvdt, false, dvdp, dvdn);
+	get_export_name(eos, specificvolume)(&T, &p, z.data(), &phase, &v.value_, v.dt_ptr, v.dp_ptr, v.dn_ptr);
+	return v;
+    }
+
+    Property zfac(double T, double p, vector1d z, int phase, bool dzdt=false, bool dzdp=false, bool dzdn=false) const {
+	activate();
+	Property Z(nc, dzdt, false, dzdp, dzdn);
+	get_export_name(eos, zfac)(&T, &p, z.data(), &phase, &Z.value_, Z.dt_ptr, Z.dp_ptr, Z.dn_ptr);
+	return Z;
+    }
+
+    VectorProperty thermo(double T, double p, vector1d z, int phase, bool dlnfugdt=false, bool dlnfugdp=false, bool dlnfugdn=false, int* ophase=nullptr, double* v=nullptr) const {
+	int metaextremum;
+	activate();
+	VectorProperty lnfug(nc, dlnfugdt, false, dlnfugdp, dlnfugdn);
+	get_export_name(eos, thermo)(&T, &p, z.data(), &phase, lnfug.value_.data(), lnfug.dt_ptr, lnfug.dp_ptr, lnfug.dn_ptr, ophase, &metaextremum, v);
+	return lnfug;
+    }
+
+    Property enthalpy(double T, double p, vector1d z, int phase, bool dhdt=false, bool dhdp=false, bool dhdn=false, int property_flag=PropertyFlag::total) const {
+	activate();
+	Property h(nc, dhdt, false, dhdp, dhdn);
+	get_export_name(eos, enthalpy)(&T, &p, z.data(), &phase, &h.value_, h.dt_ptr, h.dp_ptr, h.dn_ptr, &property_flag);
+	return h;
+    }
+
+    Property entropy(double T, double p, vector1d z, int phase, bool dsdt=false, bool dsdp=false, bool dsdn=false, int property_flag=PropertyFlag::total) const {
+	activate();
+	Property s(nc, dsdt, false, dsdp, dsdn);
+	get_export_name(eos, entropy)(&T, &p, z.data(), &phase, &s.value_, s.dt_ptr, s.dp_ptr, s.dn_ptr, &property_flag);
+	return s;
+    }
+
+    Property idealenthalpysingle(double T, int comp_idx, bool dhdt=false) const {
+	activate();
+	Property h_id(1, dhdt, false, false, false);
+	get_export_name(eos, ideal_enthalpy_single)(&T, &comp_idx, &h_id.value_, h_id.dt_ptr);
+	return h_id;
+    }
+
+    Property idealentropysingle(double T, double p, int comp_idx, bool dsdt=false, bool dsdp=false) const {
+	activate();
+	Property s_id(1, dsdt, false, dsdp, false);
+	get_export_name(eos, ideal_entropy_single)(&T, &p, &comp_idx, &s_id.value_, s_id.dt_ptr, s_id.dp_ptr);
+	return s_id;
+    }
+
+/************************************************************************/
+// ------------------------ TV-Property Interfaces -------------------- //
+/************************************************************************/
+
+    Property pressure_tv(double T, double V, vector1d n, bool dpdt=false, bool dpdv=false, bool dpdn=false,
+			int property_flag=PropertyFlag::total) const {
+	activate();
+	Property p(nc, dpdt, dpdv, false, dpdn);
+	p.value_ = get_export_name(eostv, pressure)(&T, &V, n.data(), p.dv_ptr, p.dt_ptr, nullptr, p.dn_ptr, &property_flag);
+	return p;
+    }
+
+    Property internal_energy_tv(double T, double V, vector1d n, bool dudt=false, bool dudv=false, bool dudn=false, int property_flag=PropertyFlag::total) const {
+	activate();
+	Property U(nc, dudt, dudv, false, dudn);
+	get_export_name(eostv, internal_energy_tv)(&T, &V, n.data(), &U.value_, U.dt_ptr, U.dv_ptr, U.dn_ptr, &property_flag);
+	return U;
+    }
+
+    Property entropy_tv(double T, double V, vector1d n, bool dsdt=false, bool dsdv=false, bool dsdn=false, int property_flag=PropertyFlag::total) const {
+	activate();
+	Property S(nc, dsdt, dsdv, false, dsdn);
+	get_export_name(eostv, entropy_tv)(&T, &V, n.data(), &S.value_, S.dt_ptr, S.dv_ptr, S.dn_ptr, &property_flag);
+	return S;
+    }
+
+    Property enthalpy_tv(double T, double V, vector1d n, bool dhdt=false, bool dhdv=false, bool dhdn=false, int property_flag=PropertyFlag::total) const {
+	activate();
+	Property H(nc, dhdt, dhdv, false, dhdn);
+	get_export_name(eostv, enthalpy_tv)(&T, &V, n.data(), &H.value_, H.dt_ptr, H.dv_ptr, H.dn_ptr, &property_flag);
+	return H;
+    }
+
+    Property helmholtz_tv(double T, double V, vector1d n, bool dadt=false, bool dadv=false, bool dadn=false, int property_flag=PropertyFlag::total) const {
+	activate();
+	Property A(nc, dadt, dadv, false, dadn);
+	get_export_name(eostv, free_energy_tv)(&T, &V, n.data(), &A.value_, A.dt_ptr, A.dv_ptr, A.dn_ptr, &property_flag);
+	return A;
+    }
+
+    VectorProperty chemical_potential_tv(double T, double V, vector1d n, bool dmudt=false, bool dmudv=false, bool dmudn=false,
+				    int property_flag=PropertyFlag::total) const {
+	    activate();
+	    VectorProperty mu(nc, dmudt, dmudv, false, dmudn);
+	    get_export_name(eostv, chemical_potential_tv)(&T, &V, n.data(), mu.value_.data(), mu.dt_ptr, mu.dv_ptr, mu.dn_ptr, &property_flag);
+	    return mu;
+    }
+
+    VectorProperty fugacity_tv(double T, double V, vector1d n, bool dlnphidt=false, bool dlnphidv=false, bool dlnphidn=false) const {
+	activate();
+	VectorProperty lnphi(nc, dlnphidt, dlnphidv, false, dlnphidn);
+	get_export_name(eostv, thermo_tv)(&T, &V, n.data(), lnphi.value_.data(), lnphi.dt_ptr, lnphi.dv_ptr, lnphi.dn_ptr);
+	return lnphi;
+    }
+
+/************************************************************************/
+// --------------------------- Flash interfaces ----------------------- //
+/************************************************************************/
+    FlashResult two_phase_tpflash(double T, double p, vector1d z) const {
+	activate();
+	FlashResult fr(T, p, z, "TP");
+	get_export_name(tp_solver, twophasetpflash)(&T, &p, z.data(), &fr.betaV, &fr.betaL, &fr.phase, fr.x.data(), fr.y.data());
+	return fr;
+    }
+
+    FlashResult two_phase_psflash(double p, double s, vector1d z, double T=0.) const {
+	activate();
+	FlashResult fr(T, p, z, "PS");
+	int ierr = 0;
+	get_export_name(ps_solver, twophasepsflash)(&fr.T, &p, z.data(), &fr.betaV, &fr.betaL, fr.x.data(), fr.y.data(), &s, &fr.phase, &ierr);
+	if (ierr) throw std::runtime_error("PS flash failed");
+	return fr;
+    }
+
+    FlashResult two_phase_phflash(double p, double h, vector1d z, double T=0.) const {
+	activate();
+	FlashResult fr(T, p, z, "PH");
+	int ierr = 0;
+	get_export_name(ph_solver, twophasephflash)(&fr.T, &p, z.data(), &fr.betaV, &fr.betaL, fr.x.data(), fr.y.data(), &h, &fr.phase, &ierr);
+	if (ierr) throw std::runtime_error("PH flash failed");
+	return fr;
+    }
+
+    FlashResult two_phase_uvflash(double u, double v, vector1d z, double T=0., double p=0.) const {
+	activate();
+	FlashResult fr(T, p, z, "UV");
+	get_export_name(uv_solver, twophaseuvflash)(&fr.T, &fr.p, z.data(), &fr.betaV, &fr.betaL, fr.x.data(), fr.y.data(), &u, &v, &fr.phase);
+	return fr;
+    }
+
+    FlashResult bubble_temperature(double p, vector1d z) const {
+	activate();
+	FlashResult bub = FlashResult::bub(300., p, z, "Bub_T");
+	int ierr;
+	bub.T = get_export_name(saturation, safe_bubt)(&p, z.data(), bub.y.data(), &ierr);
+	if (ierr == true_int) throw std::runtime_error("Bubble temperature calculation failed!");
+
+	return bub;
+    }
+
+    FlashResult bubble_pressure(double T, vector1d z) const {
+	activate();
+	FlashResult bub = FlashResult::bub(T, 1e5, z, "Bub_P");
+	int ierr;
+	bub.p = get_export_name(saturation, safe_bubp)(&T, z.data(), bub.y.data(), &ierr);
+	if (ierr == true_int)throw std::runtime_error("Bubble pressure calculation failed!");
+
+	return bub;
+    }
+
+    FlashResult dew_temperature(double p, vector1d z) const {
+	activate();
+	FlashResult dew = FlashResult::dew(300., p, z, "Dew_T");
+	int ierr;
+	dew.T = get_export_name(saturation, safe_dewt)(&p, dew.x.data(), z.data(), &ierr);
+	if (ierr == true_int)throw std::runtime_error("Dew temperature calculation failed!");
+
+	return dew;
+    }
+
+    FlashResult dew_pressure(double T, vector1d z) const {
+	activate();
+	FlashResult dew = FlashResult::dew(T, 1e5, z, "Dew_P");
+	int ierr;
+	dew.p = get_export_name(saturation, safe_dewp)(&T, dew.x.data(), z.data(), &ierr);
+	if (ierr == true_int)throw std::runtime_error("Dew pressure calculation failed!");
+
+	return dew;
+    }
+
+    vector1d critical(vector1d n, double T=0., double V=0., double tol=1e-7, double v_min=-1.) const {
+	activate();
+	double p = 0.;
+	vector1d tvp = {T, V, p};
+	int ierr;
+	get_export_name(critical, calccriticaltv)(&(tvp[0]), &(tvp[1]), n.data(), &ierr, &tol, (v_min > 0) ? &v_min : nullptr, &(tvp[2]));
+	if (ierr == true_int) throw std::runtime_error("Critical calculation failed!");
+
+	return tvp;
+    }
+
+    protected:
+    const int true_int;
+    size_t nc;
+    void activate() const {
+	get_export_name(thermopack_var, activate_model)(&model_idx);
+    }
+
+    private:
+    int model_idx;
+};
diff --git a/addon/cppThermopack/cppThermopack/utils.h b/addon/cppThermopack/cppThermopack/utils.h
new file mode 100644
index 00000000..1a13abe1
--- /dev/null
+++ b/addon/cppThermopack/cppThermopack/utils.h
@@ -0,0 +1,180 @@
+#pragma once
+#include <vector>
+#include <string>
+#include <iostream>
+
+using vector1d = std::vector<double>;
+using vector2d = std::vector<vector1d>;
+
+enum PropertyFlag{
+    total = 0,
+    residual,
+    ideal
+};
+
+enum Phase{
+    two = 0,
+    liq,
+    vap,
+    fake,
+    mingibbs,
+    solid
+};
+
+inline void check_has_diff(bool has_diff, std::string diff){
+    if (!has_diff) throw std::runtime_error("Property does not have differential " + diff);
+}
+
+class Property{
+    public:
+    inline operator double() const {return value_;}
+    inline double operator=(double val){value_ = val; return *this;}
+    inline double value() const {return value_;}
+    inline double dt() const {check_has_diff(dt_ptr, "dt"); return dt_;}
+    inline double dv() const {check_has_diff(dv_ptr, "dv"); return dv_;}
+    inline double dp() const {check_has_diff(dp_ptr, "dp"); return dp_;}
+    inline vector1d dn() const {check_has_diff(dn_ptr, "dn"); return dn_;}
+
+    friend std::ostream& operator<<(std::ostream& strm, const Property& self){
+	bool has_derivs = (self.dt_ptr || self.dv_ptr || self.dp_ptr || self.dn_ptr);
+	if (has_derivs) strm << "Property with value : ";
+	strm << self.value_;
+	if (has_derivs) strm << "\n";
+
+	if (self.dt_ptr) strm << "\tdt : " << self.dt_ << "\n";
+	if (self.dv_ptr) strm << "\tdv : " << self.dv_ << "\n";
+	if (self.dp_ptr) strm << "\tdp : " << self.dp_ << "\n";
+	if (self.dn_ptr){
+	    strm << "\tdn : ";
+	    for (double dni : self.dn_) strm << dni << ", ";
+	    strm << "\n";
+	}
+	return strm;
+    }
+
+    friend class Thermo;
+    private:
+    Property(size_t nc, bool dt, bool dv, bool dp, bool dn)
+	: dn_(nc, 0.)
+    {
+	dt_ptr = dt ? &dt_ : nullptr;
+	dv_ptr = dv ? &dv_ : nullptr;
+	dp_ptr = dp ? &dp_ : nullptr;
+	dn_ptr = dn ? dn_.data() : nullptr;
+    }
+
+    inline double** dt_p(){return &dt_ptr;}
+    double value_;
+    double dt_, dv_, dp_;
+    vector1d dn_;
+
+    double *dt_ptr, *dv_ptr, *dp_ptr, *dn_ptr;
+};
+
+class VectorProperty{
+    public:
+    inline operator vector1d() const {return value_;}
+    inline vector1d value() const {return value_;}
+    inline vector1d dt() const {check_has_diff(dt_ptr, "dt"); return dt_;}
+    inline vector1d dv() const {check_has_diff(dv_ptr, "dv"); return dv_;}
+    inline vector1d dp() const {check_has_diff(dp_ptr, "dp"); return dp_;}
+    inline vector2d dn() const {
+	check_has_diff(dn_ptr, "dn");
+	size_t nvals = value_.size();
+	vector2d dn_2d(nvals, vector1d(nvals, 0.));
+	for (size_t i = 0; i < nvals; i++){
+	    for (size_t j = 0; j < nvals; j++){
+		dn_2d[i][j] = dn_[j * nvals + i]; // dn_ uses Fortran-ordering internally
+	    }
+	}
+	return dn_2d;
+    }
+
+    friend std::ostream& operator<<(std::ostream& strm, const VectorProperty& self){
+	bool has_derivs = (self.dt_ptr || self.dv_ptr || self.dp_ptr || self.dn_ptr);
+	if (has_derivs) strm << "VectorProperty with value : ";
+	for (double d : self.value_) strm << d << ", ";
+	if (has_derivs) strm << "\n";
+
+	if (self.dt_ptr){strm << "\tdt : "; for (double d : self.dt_) strm << d << ", "; strm << "\n";}
+	if (self.dv_ptr){strm << "\tdv : "; for (double d : self.dv_) strm << d << ", "; strm << "\n";}
+	if (self.dp_ptr){strm << "\tdp : "; for (double d : self.dp_) strm << d << ", "; strm << "\n";}
+	if (self.dn_ptr){
+	    strm << "\tdn : [\n";
+	    for (vector1d line : self.dn()){
+		strm << "\t\t[ ";
+		for (double d : line){
+		    strm << d << ", ";
+		}
+		strm << "]\n";
+	    }
+	    strm << "\t     ]\n";
+	}
+	return strm;
+    }
+
+    friend class Thermo;
+    private:
+    VectorProperty(size_t size, bool dt, bool dv, bool dp, bool dn)
+	: value_(size, 0.), dt_(size, 0.), dv_(size, 0.), dp_(size, 0.), dn_(size * size, 0.)
+    {
+	dt_ptr = dt ? dt_.data() : nullptr;
+	dv_ptr = dv ? dv_.data() : nullptr;
+	dp_ptr = dp ? dp_.data() : nullptr;
+	dn_ptr = dn ? dn_.data() : nullptr;
+    }
+
+    vector1d value_;
+    vector1d dt_, dv_, dp_, dn_;
+    double *dt_ptr, *dv_ptr, *dp_ptr, *dn_ptr;
+
+};
+
+class FlashResult{
+    public:
+    int phase;
+    vector1d x, y, z;
+    double T, p, betaV, betaL;
+    std::string flash_type;
+
+    FlashResult(FlashResult&&) = default;
+    FlashResult(const FlashResult&) = default;
+    FlashResult& operator=(FlashResult&&) = default;
+
+    friend class Thermo;
+    friend std::ostream& operator<<(std::ostream& os, const FlashResult& self){
+	os << "FlashResult from " << self.flash_type << "-flash :\n"
+	   << "\tT : " << self.T << " K\n"
+	   << "\tp : " << self.p << " Pa\n"
+	   << "\tz : {";
+	for (const double& zi : self.z) os << zi << ", ";
+	os << "}\n"
+	   << "\tbetaL / betaV : " << self.betaL << " / " << self.betaV << "\n"
+	   << "\tx : {";
+	for (const double& xi : self.x) os << xi << ", ";
+	os << "}\n\ty : {";
+	for (const double& yi : self.y) os << yi << ", ";
+	os << "}\n";
+	return os;
+    }
+    private:
+    FlashResult(const vector1d& z_, std::string flash_type)
+	: phase{0}, x(z_.size(), 0.), y(z_.size(), 0.), z{z_}, flash_type{flash_type}
+	{}
+
+    FlashResult(double T, double p, const vector1d& z_, std::string flash_type)
+	: phase{0}, x(z_.size(), 0.), y(z_.size(), 0.), z{z_}, T{T}, p{p}, flash_type{flash_type}
+	{}
+
+    static FlashResult dew(double T, double p, const vector1d& z_, std::string spec){
+	FlashResult fr(T, p, z_, spec);
+	fr.y = z_; fr.betaV = 1.0; fr.betaL = 0.0;
+	return fr;
+    }
+
+    static FlashResult bub(double T, double p, const vector1d& z_, std::string spec){
+	FlashResult fr(T, p, z_, spec);
+	fr.x = z_; fr.betaL = 1.0; fr.betaV = 0.0;
+	return fr;
+    }
+};
diff --git a/addon/pyUtils/docs/tools.py b/addon/pyUtils/docs/tools.py
index 7d9fc66a..5e8fd3c9 100644
--- a/addon/pyUtils/docs/tools.py
+++ b/addon/pyUtils/docs/tools.py
@@ -19,7 +19,7 @@ def remove_illegal_link_chars(link_str):
     Remove characters that are illegal in url-links from the string 'link_str', and return the new string.
     """
     return (link_str.replace(' ', '-').replace('(', '').replace(')', '').replace('=', '').replace(',', '').replace(
-             "'", '').replace('.','')).replace(':', '').lower()
+             "'", '').replace('.','')).replace(':', '').replace('*', '').lower()
 
 def check_is_changed(old_file, new_file_str):
     """
diff --git a/addon/pycThermopack/.gitignore b/addon/pycThermopack/.gitignore
index 35a02010..e19f946a 100644
--- a/addon/pycThermopack/.gitignore
+++ b/addon/pycThermopack/.gitignore
@@ -3,5 +3,4 @@ LICENSE*
 build/
 thermopack/platform_specifics.py
 pyctp/platform_specifics.py
-setup.py
-pyproject.toml
\ No newline at end of file
+setup.py
\ No newline at end of file
diff --git a/addon/pycThermopack/CMakeLists.txt b/addon/pycThermopack/CMakeLists.txt
new file mode 100644
index 00000000..e69de29b
diff --git a/addon/pycThermopack/pyproject.toml b/addon/pycThermopack/pyproject.toml
new file mode 100644
index 00000000..8bb46ac7
--- /dev/null
+++ b/addon/pycThermopack/pyproject.toml
@@ -0,0 +1,40 @@
+[build-system]
+requires = ["scikit-build-core"]
+build-backend = "scikit_build_core.build"
+
+[tool.scikit-build]
+sdist.include = ["libthermopack*"]
+wheel.expand-macos-universal-tags = true
+
+[project]
+name = "thermopack"
+version = "2.3b0"
+authors = [
+  { name = "Morten Hammer", email="morten.hammer@ntnu.no" },
+]
+maintainers = [
+  { name = "Morten Hammer", email="morten.hammer@ntnu.no" },
+  { name = "Vegard Gjeldvik Jervell", email="vegard.g.jervell@ntnu.no" },
+]
+description = "Thermodynamics library for multicomponent, multiphase equilibrium calculations."
+readme = "README_pypi.md"
+requires-python = ">=3.6"
+classifiers = [
+    "Programming Language :: Python :: 3 :: Only",
+    "License :: OSI Approved :: MIT License",
+]
+keywords = ["physics", "thermodynamics", "equations_of_state", "phase_equilibria", "SAFT"]
+dependencies = ["numpy~=1.20"]
+
+[project.urls]
+"Homepage" = "https://thermotools.github.io/thermopack/"
+"Repository" = "https://github.com/thermotools/thermopack/"
+"Bug Tracker" = "https://github.com/thermotools/thermopack/issues"
+
+[project.optional-dependencies]
+test = ["pytest"]
+
+[tool.cibuildwheel]
+test-command = "pytest {project}/../pyTests"
+test-requires = ["pytest"]
+build-verbosity = 2
\ No newline at end of file
diff --git a/addon/pycThermopack/thermopack/thermo.py b/addon/pycThermopack/thermopack/thermo.py
index 95dabad4..8e5e05ae 100644
--- a/addon/pycThermopack/thermopack/thermo.py
+++ b/addon/pycThermopack/thermopack/thermo.py
@@ -1587,7 +1587,7 @@ def idealenthalpysingle(self, temp, j, dhdt=None):
         prop = utils.Property.from_return_tuple(return_tuple, (dhdt, None, None), 'tpn')
         return prop.unpack()
 
-    def idealentropysingle(self,temp,press,j,dsdt=None,dsdp=None):
+    def idealentropysingle(self, temp, press, j, dsdt=None, dsdp=None):
         """Tp-property
         Calculate specific ideal entropy
         Note that the order of the output match the default order of input for the differentials.
diff --git a/docs/_includes/sidebar.html b/docs/_includes/sidebar.html
index 4d890be4..7a54118a 100644
--- a/docs/_includes/sidebar.html
+++ b/docs/_includes/sidebar.html
@@ -1,4 +1,8 @@
-<ul>
+<!--
+Automatically generated file
+Generated by parsing : sidebar.md
+Using the file at : docs/_includes/sidebar_to_html.py
+--><ul>
 <li><a href="/thermopack/vcurrent/source_build.html">Installation</a></li>
 <li>
 	<a href="/thermopack/vcurrent/getting_started.html">Getting started</a> <button class="expand-btn">&#9668;</button>
@@ -46,6 +50,14 @@
 		</ul>
 	</ul>
 </li>
+<li>
+	<a href="/thermopack/vcurrent/getting_started_cpp.html">Getting Started - C++</a> <button class="expand-btn">&#9668;</button>
+	<ul class="submenu">
+		<ul>
+		<li><a href="/thermopack/vcurrent/getting_started_cpp.html#examples">Examples</a></li>
+		</ul>
+	</ul>
+</li>
 <li><a href="/thermopack/vcurrent/Component-name-mapping.html">Component identifiers</a></li>
 <li><a href="/thermopack/vcurrent/phase_flags.html">Phase flags</a></li>
 <li>
diff --git a/docs/_includes/sidebar.md b/docs/_includes/sidebar.md
index 1e3271f0..f81ce2c5 100644
--- a/docs/_includes/sidebar.md
+++ b/docs/_includes/sidebar.md
@@ -13,8 +13,10 @@
   * [Isolines](/thermopack/vcurrent/getting_started.html#isolines)
   * [Critical point](/thermopack/vcurrent/getting_started.html#critical-point)
 * [Advanced usage - Python](/thermopack/vcurrent/more_advanced.html)
-  * [Interaction parameters](/thermopack/vcurrent/more_advanced.html#interaction-parameters) 
-  * [Tuning cubics](/thermopack/vcurrent/more_advanced.html#tuning-cubics) 
+    * [Interaction parameters](/thermopack/vcurrent/more_advanced.html#interaction-parameters) 
+    * [Tuning cubics](/thermopack/vcurrent/more_advanced.html#tuning-cubics) 
+* [Getting Started - C++](/thermopack/vcurrent/getting_started_cpp.html)
+  * [Examples](/thermopack/vcurrent/getting_started_cpp.html#examples)
 * [Component identifiers](/thermopack/vcurrent/Component-name-mapping.html) 
 * [Phase flags](/thermopack/vcurrent/phase_flags.html)
 * [Adding new fluids](/thermopack/vcurrent/new_fluids.html)
diff --git a/docs/_includes/sidebar_2.1.0.html b/docs/_includes/sidebar_2.1.0.html
index 34280164..e3b85f68 100644
--- a/docs/_includes/sidebar_2.1.0.html
+++ b/docs/_includes/sidebar_2.1.0.html
@@ -1,4 +1,8 @@
-<ul>
+<!--
+Automatically generated file
+Generated by parsing : sidebar_2.1.0.md
+Using the file at : docs/_includes/sidebar_to_html.py
+--><ul>
 <li><a href="/thermopack/v2.1.0/installation.html">Installation</a></li>
 <li>
 	<a href="/thermopack/v2.1.0/getting_started.html">Getting started</a> <button class="expand-btn">&#9668;</button>
diff --git a/docs/_includes/sidebar_2.2.0.html b/docs/_includes/sidebar_2.2.0.html
index 002b43be..a53caa0b 100644
--- a/docs/_includes/sidebar_2.2.0.html
+++ b/docs/_includes/sidebar_2.2.0.html
@@ -1,4 +1,8 @@
-<ul>
+<!--
+Automatically generated file
+Generated by parsing : sidebar_2.2.0.md
+Using the file at : docs/_includes/sidebar_to_html.py
+--><ul>
 <li><a href="/thermopack/v2.2.0/source_build.html">Installation</a></li>
 <li>
 	<a href="/thermopack/v2.2.0/getting_started.html">Getting started</a> <button class="expand-btn">&#9668;</button>
diff --git a/docs/_includes/sidebar_to_html.py b/docs/_includes/sidebar_to_html.py
index cd81ba32..550abe1e 100644
--- a/docs/_includes/sidebar_to_html.py
+++ b/docs/_includes/sidebar_to_html.py
@@ -1,5 +1,5 @@
-import os
-import re
+import os, re
+from datetime import datetime
 
 link_pattern = r"\[([^]]+)\]\(([^)]+)\)"
 code_pattern = r'`([^`]+)`'
@@ -72,6 +72,11 @@ def write_table(file, tab, level):
             table = get_table(None, '', ifile, 0)
 
         with open(ofilename, 'w') as ofile:
+            ofile.write('<!--\n'
+                        'Automatically generated file\n'
+                        f'Generated by parsing : {ifilename}\n'
+                        'Using the file at : docs/_includes/sidebar_to_html.py\n'
+                        '-->')
             write_table(ofile, table, 0)
             ofile.write("""<script>
             // Get all expand buttons
diff --git a/docs/vCurrent/getting_started.md b/docs/vCurrent/getting_started.md
index cb0577c5..18917e03 100644
--- a/docs/vCurrent/getting_started.md
+++ b/docs/vCurrent/getting_started.md
@@ -1,7 +1,7 @@
 ---
 layout: default
 version: 
-title: Getting Started
+title: Getting Started - Python
 permalink: /vcurrent/getting_started.html
 ---
 
@@ -14,17 +14,20 @@ find the appropriate version in the sidebar on the [ThermoPack homepage.](https:
 Equations of State (EoS's) in ThermoPack are classes. To do calculations for a given mixture an EoS object must first be initialized for that mixture, as demonstrated in the [Initializing an EoS section](#initialising-an-equation-of-state). Then, a wide variety of thermodynamic computations can be done, as demonstrated in the remaining sections.
 
 ## Contents
-* [Initialising an equation of state](#initialising-an-equation-of-state)
-* [pVT properties](#pvt-properties)
-  * [Differentials](#differentials)
-* [Phase diagrams and equilibria](#phase-diagrams-and-equilibria)
-  * [Flash calculations](#flash-calculations)
-  * [Phase envelopes](#phase-envelopes)
-    * [Tp- and Tv- envelopes](#tp--and-tv--phase-envelopes)
-    * [pxy- and Txy- envelopes](#pxy--and-txy--phase-envelopes)
-  * [Dew- and bubble points](#dew--and-bubble-points)
-* [Isolines](#isolines)
-* [Critical point](#critical-point)
+- [Getting started - Python](#getting-started---python)
+  - [Contents](#contents)
+  - [Initialising an equation of state](#initialising-an-equation-of-state)
+- [Doing calculations](#doing-calculations)
+  - [pVT-properties](#pvt-properties)
+    - [Differentials](#differentials)
+  - [Phase diagrams and Equilibria](#phase-diagrams-and-equilibria)
+    - [Flash calculations](#flash-calculations)
+    - [Phase envelopes](#phase-envelopes)
+      - [Tp- and Tv- phase envelopes](#tp--and-tv--phase-envelopes)
+      - [pxy- and txy- phase envelopes](#pxy--and-txy--phase-envelopes)
+    - [Dew- and bubble points](#dew--and-bubble-points)
+  - [Isolines](#isolines)
+  - [Critical point](#critical-point)
 
 ## Initialising an equation of state
 An overview of available equations of state can be found [here](method_docs.html).
diff --git a/docs/vCurrent/getting_started_cpp.md b/docs/vCurrent/getting_started_cpp.md
new file mode 100644
index 00000000..afab70d1
--- /dev/null
+++ b/docs/vCurrent/getting_started_cpp.md
@@ -0,0 +1,85 @@
+---
+layout: default
+version: 
+title: Getting Started - C++
+permalink: /vcurrent/getting_started_cpp.html
+---
+
+The ThermoPack C++ wrapper is a header-only wrapper designed to be compatible with `stl` containers like `std::vector`, and to support use of `std::unique_ptr` and cousins. In short, it is intended to make using ThermoPack together with modern C++ as painless as possible.
+
+- [Building ThermoPack](#building-thermopack)
+- [ThermoPack headers](#thermopack-headers)
+- [Including and linking ThermoPack](#including-and-linking-thermopack)
+- [Library structure](#library-structure)
+- [Examples](#examples)
+
+# Building ThermoPack
+
+The ThermoPack C++ wrapper is a header-only wrapper, that simplifies linking to the `libthermopack` library that is built and installed with the `CMakeLists.txt` in the root directory of ThermoPack.
+
+To build the required library,
+
+```bash
+git clone https://github.com/thermotools/thermopack.git
+cd thermopack
+mkdir build
+cd build
+cmake ..
+make install
+```
+
+This should generate the file `thermopack/addon/cppThermopack/libthermopack.[so/dylib]`, which you will later be linking to.
+
+# ThermoPack headers
+
+The header files you will need are found in `thermopack/addon/cppThermoPack/cppThermoPack`. 
+
+# Including and linking ThermoPack
+
+An example `CMakeLists.txt` is found in `thermopack/addon/cppExamples`. In essence, all that is required to include ThermoPack is to add the following to your `CMakeLists.txt`:
+
+```cmake
+set(THERMOPACK_HEADERS "path/to/thermopack/addon/cppThermopack")
+set(libthermopack_path "path/to/thermopack/addon/cppThermopack/libthermopack.[so/dylib]")
+
+target_link_libraries(<my_project> ${libthermopack_path})
+target_include_directories(<my_project> ${THERMOPACK_HEADERS})
+```
+
+where `<my_project>` is some target previously defined by a call to `add_executable`, `add_library`, etc.
+
+# Library structure
+
+The ThermoPack C++ wrapper defines a class structure that is almost identical to the structure in the python wrapper. All equations of state inherit from the `Thermo` class (defined in `thermo.h`), and only implement initialization and some minor utility methods.
+
+In addition, several utility structures are defined in `utils.h`. These include the `Property` and `VectorProperty` classes, which are used to hold computed values and derivatives. The structures support implicit conversion to respectively `double` and `std::vector<double>`, such that their existence can be ignored unless you are interested in computing derivatives.
+
+The result of phase equilibrium calculations, such as flash- and saturation point calculations, are returned as a `FlashResult` object, which holds the attributes
+* `double T` - Temperature (K)
+* `double p` - Pressure (Pa)
+* `double betaV` - Vapour fraction
+* `double betaL` - Liquid fraction
+* `std::vector<double> z` - Total composition
+* `std::vector<double> x` - Liquid composition
+* `std::vector<double> y` - Vapour composition
+* `std::string flash_type` - Information about the calculation this object originated from
+
+# Examples
+
+Some sample programs using ThermoPack in C++ are found in the directory `thermopack/addon/cppExamples`. To build them, navigate to that directory, and
+
+```bash
+mkdir build
+cd build
+cmake ..
+make
+```
+
+This will build the examples as executables in the `build` directory.
+
+Currently implemented examples are
+* `basic` - Demonstrates initialization of some equations of state
+* `flashes` - Flash calculations
+* `tp_properties` - Compute properties at given temperature and pressure
+* `tv_properties` - Compute properties at given temperature and volume
+* `saturation` - Compute dew- and bubble points, and the critical point of a mixture
\ No newline at end of file
diff --git a/docs/vCurrent/thermo_methods.md b/docs/vCurrent/thermo_methods.md
index 7c61792f..ea566428 100644
--- a/docs/vCurrent/thermo_methods.md
+++ b/docs/vCurrent/thermo_methods.md
@@ -6,7 +6,7 @@ permalink: /vcurrent/thermo_methods.html
 ---
 
 <!--- 
-Generated at: 2024-06-11T11:12:19.233427
+Generated at: 2024-07-25T14:16:10.431252
 This is an auto-generated file, generated using the script at thermopack/addon/pyUtils/docs/markdown_from_docstrings.py
 The file is created by parsing the docstrings of the methods in the 
 thermo class. For instructions on how to use the parser routines, see the
@@ -98,7 +98,7 @@ The `thermo` class, found in `addon/pycThermopack/thermopack/thermo.py`, is the
     * [joule_thompson_inversion](#joule_thompson_inversionself-z-nmax1000)
   * [Utility methods](#utility-methods)
     * [acentric_factor](#acentric_factorself-i)
-    * [compmoleweight](#compmoleweightself-comp)
+    * [compmoleweight](#compmoleweightself-comp-si_unitsfalse)
     * [get_comp_name](#get_comp_nameself-index-get_comp_identifierfalse)
     * [get_comp_structure](#get_comp_structureself-comp_name)
     * [get_enthalpy_of_formation](#get_enthalpy_of_formationself-j)
@@ -112,6 +112,7 @@ The `thermo` class, found in `addon/pycThermopack/thermopack/thermo.py`, is the
     * [get_tmax](#get_tmaxself)
     * [get_tmin](#get_tminself)
     * [getcompindex](#getcompindexself-comp)
+    * [moleweight](#moleweightself-z-si_unitstrue)
     * [redefine_critical_parameters](#redefine_critical_parametersself-silenttrue-tc_initialsnone-vc_initialsnone)
     * [set_enthalpy_of_formation](#set_enthalpy_of_formationself-j-h0)
     * [set_ideal_cp](#set_ideal_cpself-j-cp_correlation_type-parameters)
@@ -124,7 +125,7 @@ The `thermo` class, found in `addon/pycThermopack/thermopack/thermo.py`, is the
   * [Internal methods](#internal-methods)
     * [\_\_del\_\_](#__del__self)
     * [\_\_init\_\_](#__init__self)
-    * [\_\_new\_\_](#__new__cls-*args-**kwargs)
+    * [\_\_new\_\_](#__new__cls-args-kwargs)
     * [_get_true_int_value](#_get_true_int_valueself)
     * [activate](#activateself)
     * [add_eos](#add_eosself)
@@ -2941,7 +2942,7 @@ Methods for setting ... and getting ...
 ### Table of contents
   * [Utility methods](#utility-methods)
     * [acentric_factor](#acentric_factorself-i)
-    * [compmoleweight](#compmoleweightself-comp)
+    * [compmoleweight](#compmoleweightself-comp-si_unitsfalse)
     * [get_comp_name](#get_comp_nameself-index-get_comp_identifierfalse)
     * [get_comp_structure](#get_comp_structureself-comp_name)
     * [get_enthalpy_of_formation](#get_enthalpy_of_formationself-j)
@@ -2955,6 +2956,7 @@ Methods for setting ... and getting ...
     * [get_tmax](#get_tmaxself)
     * [get_tmin](#get_tminself)
     * [getcompindex](#getcompindexself-comp)
+    * [moleweight](#moleweightself-z-si_unitstrue)
     * [redefine_critical_parameters](#redefine_critical_parametersself-silenttrue-tc_initialsnone-vc_initialsnone)
     * [set_enthalpy_of_formation](#set_enthalpy_of_formationself-j-h0)
     * [set_ideal_cp](#set_ideal_cpself-j-cp_correlation_type-parameters)
@@ -2971,7 +2973,7 @@ Get acentric factor of component i Args: i (int) component FORTRAN index returns
 
 &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; 
 
-### `compmoleweight(self, comp)`
+### `compmoleweight(self, comp, si_units=False)`
 Get component mole weight (g/mol)
 
 #### Args:
@@ -2980,6 +2982,10 @@ Get component mole weight (g/mol)
 
 &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;  Component FORTRAN index
 
+&nbsp;&nbsp;&nbsp;&nbsp; **si_units (bool, optional) :** 
+
+&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;  true for output in kg/mol, false for g/mol
+
 &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; 
 
 #### Returns:
@@ -3191,6 +3197,29 @@ Get component index
 
 &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; 
 
+### `moleweight(self, z, si_units=True)`
+Get mole weight (kg/mol)
+
+#### Args:
+
+&nbsp;&nbsp;&nbsp;&nbsp; **z (array_like):** 
+
+&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;  Molar composition
+
+&nbsp;&nbsp;&nbsp;&nbsp; **si_units (bool, optional) :** 
+
+&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;  true for output in kg/mol, false for g/mol
+
+&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; 
+
+#### Returns:
+
+&nbsp;&nbsp;&nbsp;&nbsp; **float:** 
+
+&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;  mixture mole weight (kg/mol)
+
+&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; 
+
 ### `redefine_critical_parameters(self, silent=True, Tc_initials=None, vc_initials=None)`
 Recalculate critical properties of pure fluids
 
@@ -3326,7 +3355,7 @@ Methods for handling communication with the Fortran library.
   * [Internal methods](#internal-methods)
     * [\_\_del\_\_](#__del__self)
     * [\_\_init\_\_](#__init__self)
-    * [\_\_new\_\_](#__new__cls-*args-**kwargs)
+    * [\_\_new\_\_](#__new__cls-args-kwargs)
     * [_get_true_int_value](#_get_true_int_valueself)
     * [activate](#activateself)
     * [add_eos](#add_eosself)
diff --git a/external/lapack b/external/lapack
new file mode 160000
index 00000000..3c351aad
--- /dev/null
+++ b/external/lapack
@@ -0,0 +1,1 @@
+Subproject commit 3c351aad417ae1ad8552bc742a6f039d6a298f2c
diff --git a/external/pFUnit b/external/pFUnit
new file mode 160000
index 00000000..26dadb11
--- /dev/null
+++ b/external/pFUnit
@@ -0,0 +1,1 @@
+Subproject commit 26dadb1157819ea1bd9c355c60ed52f42dd36432
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 00000000..2f1dc662
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,42 @@
+file(GLOB SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/*.f90)
+list(REMOVE_ITEM SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/thermopack.f90)
+set(TREND_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../addon/trend_interface/include)
+set(TREND ${TREND_DIR}/trend_interface.f95)
+
+if(NOT MSVC)
+    find_package(LAPACK REQUIRED)
+endif()
+
+# Define the executable in terms of the source files
+add_library(libthermopack STATIC ${SOURCES})
+if(MSVC)
+    add_library(thermopack SHARED ${SOURCES})
+    set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} /libs:static /threads")
+    set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} /DEF:${CMAKE_SOURCE_DIR}/MSVStudio/thermopack.def")
+else(MSVC)
+    if(APPLE)
+        set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -exported_symbols_list ${CMAKE_SOURCE_DIR}/libthermopack_export.symbols")
+    else(APPLE)
+        set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -Wl,--version-script=${CMAKE_SOURCE_DIR}/libthermopack_export.version")
+    endif(APPLE)
+    add_library(thermopack SHARED $<TARGET_OBJECTS:libthermopack>)
+endif(MSVC)
+
+set_target_properties(libthermopack PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+target_include_directories(libthermopack PUBLIC ${TREND_DIR})
+set_target_properties(thermopack PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+target_include_directories(thermopack PUBLIC ${TREND_DIR})
+
+install(TARGETS thermopack DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/../addon/pycThermopack/thermopack)
+install(TARGETS thermopack DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/../addon/cppThermopack)
+
+add_executable(run_thermopack ${CMAKE_CURRENT_SOURCE_DIR}/thermopack.f90)
+target_include_directories(run_thermopack PUBLIC ${CMAKE_CURRENT_BINARY_DIR})
+
+if(NOT MSVC)
+    target_link_libraries(thermopack ${LAPACK_LIBRARIES})
+    target_link_libraries(run_thermopack libthermopack ${LAPACK_LIBRARIES})
+else()
+    target_link_libraries(thermopack PRIVATE lapack blas)
+    target_link_libraries(run_thermopack libthermopack lapack blas)
+endif()
diff --git a/thermopack-config.cmake b/thermopack-config.cmake
new file mode 100644
index 00000000..11b3a081
--- /dev/null
+++ b/thermopack-config.cmake
@@ -0,0 +1,36 @@
+# Check if thermopack has been compiled and installed to the proper location for inclusion in other projects.
+# Result variables:
+#   THERMOPACK_INSTALLED : "TRUE" if thermopack is installed, "FALSE" otherwise
+#   THERMOPACK_ROOT : Path to directory containing this file
+#   THERMOPACK_LIB : Path to thermopack dynamic library
+#   THERMOPACK_INCLUDE : Path to thermopack C++ headers
+#   thermopack : Imported shared library target with headers
+
+
+message(STATUS "ThermoPack found at: ${CMAKE_CURRENT_LIST_DIR}")
+set(THERMOPACK_LIB ${CMAKE_CURRENT_LIST_DIR}/addon/cppThermopack/libthermopack${CMAKE_SHARED_LIBRARY_SUFFIX})
+set(THERMOPACK_ROOT ${CMAKE_CURRENT_LIST_DIR})
+
+if(NOT EXISTS ${THERMOPACK_LIB})
+    unset(THERMOPACK_LIB)
+    set(THERMOPACK_INSTALLED FALSE)
+    message(STATUS "ThermoPack has not been properly installed.")
+    return()
+endif()
+
+set(THERMOPACK_FOUND TRUE)
+set(THERMOPACK_INSTALLED TRUE)
+set(THERMOPACK_INCLUDE ${CMAKE_CURRENT_LIST_DIR}/addon/cppThermopack)
+
+message(STATUS "Setting ThermoPack cmake variables ...")
+message(STATUS "THERMOPACK_ROOT: ${THERMOPACK_ROOT}")
+message(STATUS "THERMOPACK_LIB: ${THERMOPACK_LIB}")
+message(STATUS "THERMOPACK_INCLUDE: ${THERMOPACK_INCLUDE}")
+
+if(NOT TARGET thermopack)
+    add_library(thermopack SHARED IMPORTED)
+    set_target_properties(thermopack PROPERTIES 
+                            IMPORTED_LOCATION ${THERMOPACK_LIB}
+                            INTERFACE_INCLUDE_DIRECTORIES ${THERMOPACK_INCLUDE})
+    message(STATUS "ThermoPack exported target: thermopack")
+endif()
\ No newline at end of file


Test patch:
diff --git a/unittests/CMakeLists.txt b/unittests/CMakeLists.txt
new file mode 100644
index 00000000..e54361ae
--- /dev/null
+++ b/unittests/CMakeLists.txt
@@ -0,0 +1,8 @@
+
+set(CMAKE_Fortran_FLAGS ${tp_fortran_flags})
+
+enable_testing()
+file(GLOB UNITTEST_SRC ${CMAKE_CURRENT_SOURCE_DIR}/*.pf)
+list(REMOVE_ITEM UNITTEST_SRC ${CMAKE_CURRENT_SOURCE_DIR}/test_trend.pf)
+add_pfunit_ctest(unittests TEST_SOURCES ${UNITTEST_SRC} LINK_LIBRARIES thermopack)
+target_include_directories(unittests PRIVATE ${CMAKE_BINARY_DIR}/thermopack)
\ No newline at end of file
diff --git a/unittests/test_critical.pf b/unittests/test_critical.pf
index e8e379c1..5bb52655 100644
--- a/unittests/test_critical.pf
+++ b/unittests/test_critical.pf
@@ -67,7 +67,7 @@ contains
 
     T = -1.0
     v = -1.0
-    tol = 1.0e-7
+    tol = 1.0e-8
     call calcCriticalTV(t,v,Z,ierr,tol)
     @assertTrue(ierr == 0)
     p = pressure(T,v,Z)