Repository: thermotools/thermopack (language: python)

Issue:
Initial prototype / proof-of-concept for PyBind11/Fortwrap wrapping
The current `pybind11_bindings.cpp` contains minor manual tweaks that are necessary and have not yet been automated. Specifically: The lambdas exposed to python take `std::vector`s as arguments (such that python can pass a `list`), and forwards `std::vector.data()` to the internal C++ method. This is necessary for non-modifiable, obligatory arguments that are arrays.

As a side note: Would it maybe be cleaner for the repo structure to have the branch `beta_v3/main`, such that branches related to `beta_v3` can be placed under e.g. `beta_v3/my_stuff`? I think that would make it easier to see what branches are related to the current version, and which are related to v3.


Discussion / hints:

Fix patch:
diff --git a/prototype/.gitignore b/prototype/.gitignore
new file mode 100644
index 00000000..ac8e3b43
--- /dev/null
+++ b/prototype/.gitignore
@@ -0,0 +1,10 @@
+# Ignore binaries, mod files and cmake/make output
+*.mod
+*.o
+*.a
+*.so
+build/
+
+# Ignore executables
+cpp_main
+fortran_main
\ No newline at end of file
diff --git a/prototype/CMakeLists.txt b/prototype/CMakeLists.txt
new file mode 100644
index 00000000..07ec6530
--- /dev/null
+++ b/prototype/CMakeLists.txt
@@ -0,0 +1,29 @@
+cmake_minimum_required(VERSION 3.12)
+cmake_policy(SET CMP0127 NEW)
+cmake_policy(SET CMP0048 NEW)
+
+project(DemoPack LANGUAGES Fortran CXX)
+
+if(APPLE)
+  set(PYBIND11_ROOT_DIR /usr/local/include/pybind11)
+  set(CMAKE_CXX_STANDARD 17)
+  set(CXX "/opt/homebrew/bin/g++-13")
+  set(CMAKE_FORTRAN_COMPILER /opt/homebrew/bin/gfortran)
+  # enable_language(CXX)
+  # set(CMAKE_CXX_COMPILER /opt/homebrew/bin/g++-13)
+  set(CMAKE_OSX_ARCHITECTURES "arm64" CACHE STRING "The OSX architecture")
+  set(CMAKE_CXX_FLAGS "-fPIC -std=c++17 -arch arm64 -Wfatal-errors -Wno-unused-command-line-argument -Wl,-no_compact_unwind")
+else()
+  set(PYBIND11_ROOT_DIR ${PROJECT_SOURCE_DIR}/pybind11)
+  set(CMAKE_FORTRAN_COMPILER gfortran)
+  set(CMAKE_CXX_FLAGS "-fPIC -std=c++17 -Wfatal-errors")
+  set(CMAKE_Fortran_FLAGS "-fPIC")
+endif(APPLE)
+
+add_subdirectory(./src fortran_lib)
+add_subdirectory(${PYBIND11_ROOT_DIR} pybind11)
+add_subdirectory(./cpp_wrapper cpp_wrapper)
+add_subdirectory(./fortran_test fortran_test)
+add_subdirectory(./cpp_test cpp_test)
+file(COPY py_wrapper DESTINATION "${CMAKE_CURRENT_BINARY_DIR}")
+file(COPY pytest.py DESTINATION "${CMAKE_CURRENT_BINARY_DIR}")
diff --git a/prototype/FortWrapOptions.txt b/prototype/FortWrapOptions.txt
new file mode 100644
index 00000000..c0d6d0fc
--- /dev/null
+++ b/prototype/FortWrapOptions.txt
@@ -0,0 +1,6 @@
+--no-vector
+-g
+-d ../cpp_wrapper
+-pybind11 libdemo
+
+%exclude_files CompData
\ No newline at end of file
diff --git a/prototype/README.md b/prototype/README.md
new file mode 100644
index 00000000..b822e585
--- /dev/null
+++ b/prototype/README.md
@@ -0,0 +1,152 @@
+# ThermoPack v3.0.0 Prototyping
+
+Here, a proof of concept is being developed to ensure that fortwrap is capable of wrapping the functionality we are looking
+to implement in ThermoPack v3, specifically the class structure.
+
+FortWrap is a pure python parser that generates a lightweight C++ wrapper, with pybind11 bindings for Fortran source code. 
+To install it: 
+
+```
+git clone https://github.com/thermotools/fortwrap.git
+git switch -c iso_c_bindings # The current working branch
+pip install .
+```
+
+The test case is currently:
+
+* `constants.f90` : Module that is `use`'d by all other modules, but is not wrapped.
+* `BaseEoS` : Abstract class 
+   * Inheriting classes must implement `Fres`
+   * Implements `Fideal`, dummy function for ideal gas residual helmholtz energy
+   * Implements `pressure`, which calls `Fres` and `Fideal`
+* `NotIdGas` : In `directeos.f90` - Concrete class inheriting from `BaseEoS`
+* `VariantEoS` : In `Variants.f90` - Abstract class inheriting from `BaseEoS`
+   * Inheriting classes must override `internal_comp`
+   * Also contains two parameters not found in `BaseEoS` (`param1` and `param2`)
+   * Implements `variant_common_comp`
+   * `Variant1` : Concrete class inheriting from `VariantEoS`
+     * Implements `internal_comp`
+     * Implements `Fres`, which calls `internal_comp` and `variant_common_comp`
+   * `Variant2` : Concrete class inheriting from `VariantEoS`
+     * Implements `internal_comp`, but differently than `Variant1`
+     * Implements `Fres`, which calls `internal_comp` and `variant_common_comp`
+* A (very simple/minimal) compdatadb, with utils, such that EoS objects can be initialized from identifier strings
+  * A python script in `fluids` is used to manipulate the database, but currently does not generate the fortran-side database
+  * `Variant1` and `Variant2` have been tested to initialize correctly for different species and parameter references
+  * Wrapping of these constructors has not been attempted yet.
+
+The wrapping procedure (e.g. running `build.sh`)
+ * Generates `ISO_C_BINDINGS` wrappers for the fortran methods
+ * Generates a C++ wrapper, with pybind11 bindings. Wrapper code is placed in `prototype/cpp_wrapper`
+ * The Python wrapper must currently be written manually, but the api is equivalent to the current thermopack api
+   * The python wrapper is found in `prototype/py_wrapper`
+
+# Status
+
+Current status:
+ * Pybind11 bindings must use `std::vector<T>` for list input from python, while the C++ wrapper uses arrays for input variables. Therefore, the lambdas must pass `var.data()` to the internal C++ method. This has been tested manually, but not automated yet. 
+ * See above for status on the `compdatadb` and initializing from identifier strings
+ * Something weird is going on with the top-level `CMakeLists.txt`, which causes it to constantly think it has changed the `CMAKE_CXX_COMPILER` and rerun indefinitely.
+   * To remedy, run `cmake .. -D CMAKE_CXX_COMPILER=/opt/homebrew/bin/g++-13`
+ * Also, extra compiler flags are being added somewhere, causing the compiler to use `-std=c++11`, not `-std=c++17`, and attempt to compile for `x86` instead of `arm64`.
+   * This is also fixed by explicitly specifying the compiler when calling `cmake` (the above fix).
+
+Moving forward:
+
+ * Implement Fortran-side initialisation from a struct that holds parameters.
+   * Partially complete - 
+     * Wrapping of these constructors has not been tested.
+     * Having several constructors (C++ side polymorphic constructors) has not been attempted yet.
+ * Implement internal computations that are not wrapped, to ensure that it is possible to make a minimal wrapper.
+ * Implement an external generic method (e.g. numerical solver) that accepts a function pointer of some sort, that can be used on a type bound procedure.
+   * For example: A method that takes `eos.pressure` as input, and solves `p(T, V, n) = p0` for e.g. `V`.
+   * This method does not need to be wrapped.
+
+To play around:
+ * `cpp_test/main.cpp`, `fortran_test/main.f90`, and `pytest.py` contain working examples of initialising various models and running dummy-computations in each layer.
+ * To build the project, navigate to the `prototype` directory (this directory), and
+
+```
+mkdir build
+cmake .. && make
+```
+
+  * This will build the files
+    * `build/fortran_lib/libdemopack.a` - The Fortran library
+    * `build/cpp_wrapper/libdemopack_cpp.a` - The C++ library (without python bindings)
+    * `build/cpp_wrapper/libdemo-cpython-<python_version>-<system>.so`
+  * The build system makes some assumptions regarding where to find `pybind11` (set in `PYBIND11_ROOT_DIR` of `prototype/CMakeLists.txt`) as well as the C++ compiler (also set in `prototype/CMakeLists.txt`). 
+    * It should not be necessary to modify the `CMakeLists.txt` files in the `prototype/src` or `prototype/cpp_wrapper` directories.
+  * To `import` the library to python, the import path must point to the `libdemo-cpython-<python_version>-<system>.so` file, and the imported module is called `libdemo`, such as
+
+```
+from build.cpp_wrapper import libdemo #Importing from the `prototype` directory
+```
+
+# Interfacing
+
+There are essentially two interfaces that are handled here `Fortran <-> C++` and `C++ <-> Python`. These interfaces both have
+some requirements that must be fulfilled, which results in some translation being necessary within the C++ wrapper. Currently,
+an attempt is being made to ensure that all this translation happens within lambda functions in the pybind11 bindings.
+
+The result is that the C++ wrapper is closer to Fortran than python, and that the C++ wrapper uses raw pointers
+and arrays, rather than stl containers.
+
+In short: Fortran has no concept of stl containers, and works only with pointers. Python (through pybind11) on the other hand works very well with
+e.g. `std::vector`, and has no concept of pointers.
+
+Parameters that are fortran `optional`, must be `pointer to pointer` on the C++ side.
+If the parameter is "not" passed, the innermost of the pointers is set to `nullptr`, such that fortran recieves a pointer to 
+`nullptr`, and behaves as if the parameter was not recieved :
+
+```
+// In C++
+float* val = nullptr;
+float** val_p = &val;
+
+! In Fortran
+present(val_p) ! Returns .false.
+```
+
+Parameters that are fortran `intent(inout)` or `intent(out)` must be `pointer` on the C++ side, and `np.array` on the python side.
+Using pybind11, a pointer to the underlying memory of the `np.array` is retrieved and passed to the C++ wrapper, which forwards 
+the pointer to Fortran.
+
+| Modifiable | Optional | Python                   | C++                        | Fortran  |
+|------------|----------|--------------------------|----------------------------|----------|
+| No         | No       | float                    | float                      | real     |
+| No         | No       | list                     | std::vector => array       | array    |
+| ---        | ---      | ------                   | ------                     | -----    |
+| Yes        | No       | float => np.array<float> | np.array<float> => float*  | real     |
+| Yes        | No       | list[T] => np.array<T>   | np.array<T> => T*          | real     |
+| ---        | ---      | ------                   | ------                     | -----    |
+| Yes        | Yes      | float => np.array<float> | np.array<float> => float** | real     |
+| Yes        | Yes      | list[T] => np.array<T>   | np.array<T> => T**         | real     |
+
+# Current workflow
+
+First: Navigate to the top-level directory of the `fortwrap` project (where `setup.py` is found) and `pip install -e .`
+
+Then: Navigate to the `thermopack/prototype` directory (this directory, from now referred to as `demo`). 
+
+Run the script `build.sh`, which does the following:
+
+ * Compile the Fortran source
+ * Build the archive `libdemo_fortran.a`
+ * Move the archive and `.mod` files to the `demo/wrappers` directory
+ * Run `python -m fortwrap -g -i ../FortWrapOptions.txt` from the `demo/src` directory
+ * Compile the resulting `FortranISOWrappers.f90` and add it to the `libdemo_fortran.a` archive
+ * Compile the `.cpp` files in the `demo/wrappers` directory
+ * Compile `demo/main.cpp`
+
+To compile the C++ wrappers with pybind11, you must additionally run `cmake .. && make` from the `demo/wrappers/build` directory.
+
+The `demo/wrappers/CMakeLists.txt` does the following:
+
+ * Define the `pybind_module` target
+ * Add the `demo/wrappers/.cpp` files to the target
+    * It is the binding module `pybind11_bindings.cpp` that is used, not `bindings.cpp`
+ * Link to the archive `libdemo_fortran.a`
+ * Link to `libgfortran.dylib`
+
+The scripts `rebuild_and_run.sh` and `rebuild_fortran_prog.sh` will build (and run) the files `main.cpp` and `main.f90`.
\ No newline at end of file
diff --git a/prototype/cpp_wrapper/BaseEos.cpp b/prototype/cpp_wrapper/BaseEos.cpp
new file mode 100644
index 00000000..ef6d6fd2
--- /dev/null
+++ b/prototype/cpp_wrapper/BaseEos.cpp
@@ -0,0 +1,28 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#include <cstring> // For strcpy
+#include "BaseEos.h"
+
+void BaseEos::BaseEos_init(const char* ident) {
+  // Create C array for Fortran input string data
+  char ident_c__[10+1];
+  if (ident) {
+    strncpy(ident_c__, ident, 10+1); ident_c__[10] = 0; // strncpy protects in case ident is too long
+    for (size_t i=strlen(ident_c__); i<10+1; i++) ident_c__[i] = ' '; // Add whitespace for Fortran
+  }
+  base_eos__baseeos_init_wrap(data_ptr, ident ? ident_c__ : nullptr);
+}
+
+void BaseEos::Fideal(float T, float V, const float n[], float* Fid, float** Ft, float** Fv, float* Fn[]) {
+  base_eos__fideal_wrap(data_ptr, &T, &V, n, Fid, Ft, Fv, Fn);
+}
+
+void BaseEos::pressure(float T, float V, const float n[], float* p) {
+  base_eos__pressure_wrap(data_ptr, &T, &V, n, p);
+}
+
+void BaseEos::set_Tc(float Tc) {
+  base_eos__set_tc_wrap(data_ptr, &Tc);
+}
+
diff --git a/prototype/cpp_wrapper/BaseEos.h b/prototype/cpp_wrapper/BaseEos.h
new file mode 100644
index 00000000..d6a6a9f6
--- /dev/null
+++ b/prototype/cpp_wrapper/BaseEos.h
@@ -0,0 +1,51 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#ifndef BASEEOS_H_
+#define BASEEOS_H_
+
+#ifdef _MSC_VER
+#define DLLEXPORT __declspec(dllexport)
+#else
+#define DLLEXPORT
+#endif
+
+
+#ifndef SWIG // Protect declarations from SWIG
+#include <cstdlib>
+#include "InterfaceDefs.h"
+
+extern "C" {
+  void base_eos__baseeos_init_wrap(ADDRESS c_this, const char* ident);
+  void base_eos__fideal_wrap(ADDRESS c_this, float* T, float* V, const float n[], float* Fid, float** Ft, float** Fv, float* Fn[]);
+  void base_eos__pressure_wrap(ADDRESS c_this, float* T, float* V, const float n[], float* p);
+  void base_eos__set_tc_wrap(ADDRESS c_this, float* Tc);
+}
+#endif // SWIG
+
+/**
+ * \brief Base EoS class
+ */
+class DLLEXPORT BaseEos {
+
+protected:
+  // BaseEos can not be instantiated
+  BaseEos() {}
+
+public:
+  virtual ~BaseEos() {}
+
+  void BaseEos_init(const char* ident);
+
+  void Fideal(float T, float V, const float n[], float* Fid, float** Ft=nullptr, float** Fv=nullptr, float* Fn[]=nullptr);
+
+  void pressure(float T, float V, const float n[], float* p);
+
+  void set_Tc(float Tc);
+
+  virtual void Fres(float T, float V, const float n[], float* Fres, float** Ft=nullptr, float** Fv=nullptr, float* Fn[]=nullptr) = 0;
+
+  ADDRESS data_ptr;
+};
+
+#endif /* BASEEOS_H_ */
diff --git a/prototype/cpp_wrapper/CMakeLists.txt b/prototype/cpp_wrapper/CMakeLists.txt
new file mode 100644
index 00000000..0a8e40e2
--- /dev/null
+++ b/prototype/cpp_wrapper/CMakeLists.txt
@@ -0,0 +1,14 @@
+file(GLOB cpp_wrapper_src *.cpp *.f90)
+list(REMOVE_ITEM cpp_wrapper_src ${CMAKE_CURRENT_SOURCE_DIR}/pybind11_bindings.cpp)
+
+add_library(demopack_cpp ${cpp_wrapper_src})
+target_include_directories(demopack_cpp PUBLIC ../build/fortran_lib)
+target_link_libraries(demopack_cpp PUBLIC demopack)
+
+pybind11_add_module(libdemo pybind11_bindings.cpp)
+target_link_libraries(libdemo PUBLIC demopack_cpp demopack)
+
+add_custom_command(
+  TARGET libdemo POST_BUILD
+  COMMAND ${CMAKE_COMMAND} -E copy "${CMAKE_CURRENT_BINARY_DIR}/libdemo.*.so"
+  "${CMAKE_CURRENT_BINARY_DIR}/../py_wrapper/")
diff --git a/prototype/cpp_wrapper/FortWrap.h b/prototype/cpp_wrapper/FortWrap.h
new file mode 100644
index 00000000..0f355e8e
--- /dev/null
+++ b/prototype/cpp_wrapper/FortWrap.h
@@ -0,0 +1,9 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#include "VariantEoS.h"
+#include "Variant1.h"
+#include "Variant2.h"
+#include "NotIdGas.h"
+#include "BaseEos.h"
+
diff --git a/prototype/cpp_wrapper/FortranISOWrappers.f90 b/prototype/cpp_wrapper/FortranISOWrappers.f90
new file mode 100644
index 00000000..1d4e5b04
--- /dev/null
+++ b/prototype/cpp_wrapper/FortranISOWrappers.f90
@@ -0,0 +1,336 @@
+MODULE constants_wrap_
+
+  USE ISO_C_BINDING
+  IMPLICIT NONE
+
+CONTAINS
+
+  ! C binding wrappers:
+
+END MODULE constants_wrap_
+
+
+MODULE base_eos_wrap_
+
+  USE ISO_C_BINDING
+  USE base_eos
+  IMPLICIT NONE
+
+  TYPE BaseEos_container_
+    CLASS (BaseEos), ALLOCATABLE :: c
+  END TYPE BaseEos_container_
+
+CONTAINS
+
+  ! C binding wrappers:
+
+  SUBROUTINE base_eos__baseeos_init_wrap(c_this, ident) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    TYPE(C_PTR), VALUE :: ident
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CHARACTER(10), POINTER :: ident__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    CALL C_F_POINTER(ident, ident__p)
+    CALL BaseEos_init(c_this__p%c, ident__p)
+  END SUBROUTINE base_eos__baseeos_init_wrap
+
+  SUBROUTINE base_eos__fideal_wrap(c_this, T, V, n, Fid, Ft, Fv, Fn) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: T
+    REAL(C_FLOAT) :: V
+    REAL(C_FLOAT), DIMENSION(*) :: n
+    REAL(C_FLOAT) :: Fid
+    REAL(C_FLOAT), OPTIONAL :: Ft
+    REAL(C_FLOAT), OPTIONAL :: Fv
+    REAL(C_FLOAT), OPTIONAL, DIMENSION(*) :: Fn
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    CALL c_this__p%c%Fideal( T, V, n, Fid, Ft, Fv, Fn)
+  END SUBROUTINE base_eos__fideal_wrap
+
+  SUBROUTINE base_eos__pressure_wrap(c_this, T, V, n, p) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: T
+    REAL(C_FLOAT) :: V
+    REAL(C_FLOAT), DIMENSION(*) :: n
+    REAL(C_FLOAT) :: p
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    CALL c_this__p%c%pressure( T, V, n, p)
+  END SUBROUTINE base_eos__pressure_wrap
+
+  SUBROUTINE base_eos__set_tc_wrap(c_this, Tc) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: Tc
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    CALL c_this__p%c%set_Tc( Tc)
+  END SUBROUTINE base_eos__set_tc_wrap
+
+END MODULE base_eos_wrap_
+
+
+MODULE notidgas_mod_wrap_
+
+  USE ISO_C_BINDING
+  USE notidgas_mod
+  USE base_eos_wrap_
+  USE base_eos
+  IMPLICIT NONE
+
+CONTAINS
+
+  ! Derived type allocate and deallocate functions:
+
+  SUBROUTINE allocate_notidgas_(NotIdGas_cptr) BIND(C)
+    TYPE (C_PTR) :: NotIdGas_cptr
+
+    TYPE (BaseEos_container_), POINTER :: NotIdGas_fptr
+
+    ALLOCATE( NotIdGas_fptr )
+    ALLOCATE( NotIdGas::NotIdGas_fptr%c )
+    NotIdGas_cptr = C_LOC(NotIdGas_fptr)
+  END SUBROUTINE allocate_notidgas_
+
+  SUBROUTINE deallocate_notidgas_(NotIdGas_cptr) BIND(C)
+    TYPE (C_PTR), VALUE :: NotIdGas_cptr
+
+    TYPE (BaseEos_container_), POINTER :: NotIdGas_fptr
+
+    CALL C_F_POINTER(NotIdGas_cptr, NotIdGas_fptr)
+    DEALLOCATE( NotIdGas_fptr )
+  END SUBROUTINE deallocate_notidgas_
+
+  ! C binding wrappers:
+
+  SUBROUTINE notidgas_mod__notidgas_ctor_sub_wrap(instance, ident, nc, Tc, V&
+    &c) BIND(C)
+    TYPE(C_PTR), VALUE :: instance
+    INTEGER(C_INT) :: ident
+    INTEGER(C_INT) :: nc
+    REAL(C_FLOAT) :: Tc
+    REAL(C_FLOAT) :: Vc
+    TYPE(BaseEos_container_), POINTER :: instance__p
+    CALL C_F_POINTER(instance, instance__p)
+    SELECT TYPE (instance__pp => instance__p%c)
+    CLASS IS (NotIdGas)
+      CALL NotIdGas_ctor_sub(instance__pp, ident, nc, Tc, Vc)
+    END SELECT
+  END SUBROUTINE notidgas_mod__notidgas_ctor_sub_wrap
+
+  SUBROUTINE notidgas_mod__notidgas_fres_wrap(c_this, T, V, n, Fres, Ft, Fv,&
+    & Fn) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: T
+    REAL(C_FLOAT) :: V
+    REAL(C_FLOAT), DIMENSION(*) :: n
+    REAL(C_FLOAT) :: Fres
+    REAL(C_FLOAT), OPTIONAL :: Ft
+    REAL(C_FLOAT), OPTIONAL :: Fv
+    REAL(C_FLOAT), OPTIONAL, DIMENSION(*) :: Fn
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (NotIdGas)
+      CALL c_this__pp%Fres( T, V, n, Fres, Ft, Fv, Fn)
+    END SELECT
+  END SUBROUTINE notidgas_mod__notidgas_fres_wrap
+
+END MODULE notidgas_mod_wrap_
+
+
+MODULE variants_wrap_
+
+  USE ISO_C_BINDING
+  USE base_eos_wrap_
+  USE base_eos
+  USE variants
+  IMPLICIT NONE
+
+CONTAINS
+
+  ! Derived type allocate and deallocate functions:
+
+  SUBROUTINE allocate_variant1_(Variant1_cptr) BIND(C)
+    TYPE (C_PTR) :: Variant1_cptr
+
+    TYPE (BaseEos_container_), POINTER :: Variant1_fptr
+
+    ALLOCATE( Variant1_fptr )
+    ALLOCATE( Variant1::Variant1_fptr%c )
+    Variant1_cptr = C_LOC(Variant1_fptr)
+  END SUBROUTINE allocate_variant1_
+
+  SUBROUTINE deallocate_variant1_(Variant1_cptr) BIND(C)
+    TYPE (C_PTR), VALUE :: Variant1_cptr
+
+    TYPE (BaseEos_container_), POINTER :: Variant1_fptr
+
+    CALL C_F_POINTER(Variant1_cptr, Variant1_fptr)
+    DEALLOCATE( Variant1_fptr )
+  END SUBROUTINE deallocate_variant1_
+
+  ! Derived type allocate and deallocate functions:
+
+  SUBROUTINE allocate_variant2_(Variant2_cptr) BIND(C)
+    TYPE (C_PTR) :: Variant2_cptr
+
+    TYPE (BaseEos_container_), POINTER :: Variant2_fptr
+
+    ALLOCATE( Variant2_fptr )
+    ALLOCATE( Variant2::Variant2_fptr%c )
+    Variant2_cptr = C_LOC(Variant2_fptr)
+  END SUBROUTINE allocate_variant2_
+
+  SUBROUTINE deallocate_variant2_(Variant2_cptr) BIND(C)
+    TYPE (C_PTR), VALUE :: Variant2_cptr
+
+    TYPE (BaseEos_container_), POINTER :: Variant2_fptr
+
+    CALL C_F_POINTER(Variant2_cptr, Variant2_fptr)
+    DEALLOCATE( Variant2_fptr )
+  END SUBROUTINE deallocate_variant2_
+
+  ! C binding wrappers:
+
+  SUBROUTINE variants__variant_common_comp_wrap(c_this, mixed_common) BIND(C&
+    &)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: mixed_common
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (VariantEoS)
+      CALL c_this__pp%variant_common_comp( mixed_common)
+    END SELECT
+  END SUBROUTINE variants__variant_common_comp_wrap
+
+  SUBROUTINE variants__variant1_ctor_sub_wrap(c_this, ident, nc, Tc, Vc, var&
+    &1, var2) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    INTEGER(C_INT) :: ident
+    INTEGER(C_INT) :: nc
+    REAL(C_FLOAT) :: Tc
+    REAL(C_FLOAT) :: Vc
+    INTEGER(C_INT) :: var1
+    REAL(C_FLOAT) :: var2
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (Variant1)
+      CALL Variant1_ctor_sub(c_this__pp, ident, nc, Tc, Vc, var1, var2)
+    END SELECT
+  END SUBROUTINE variants__variant1_ctor_sub_wrap
+
+  SUBROUTINE variants__variant1_db_ctor_sub_wrap(c_this, ident) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    TYPE(C_PTR), VALUE :: ident
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CHARACTER(10), POINTER :: ident__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    CALL C_F_POINTER(ident, ident__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (Variant1)
+      CALL Variant1_db_ctor_sub(c_this__pp, ident=ident__p)
+    END SELECT
+  END SUBROUTINE variants__variant1_db_ctor_sub_wrap
+
+  SUBROUTINE variants__variant1_internal_comp_wrap(c_this, T, V, n, computed&
+    &) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: T
+    REAL(C_FLOAT) :: V
+    REAL(C_FLOAT), DIMENSION(*) :: n
+    REAL(C_FLOAT) :: computed
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (Variant1)
+      CALL c_this__pp%internal_comp( T, V, n, computed)
+    END SELECT
+  END SUBROUTINE variants__variant1_internal_comp_wrap
+
+  SUBROUTINE variants__variant1_fres_wrap(c_this, T, V, n, Fres, Ft, Fv, Fn)&
+    & BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: T
+    REAL(C_FLOAT) :: V
+    REAL(C_FLOAT), DIMENSION(*) :: n
+    REAL(C_FLOAT) :: Fres
+    REAL(C_FLOAT), OPTIONAL :: Ft
+    REAL(C_FLOAT), OPTIONAL :: Fv
+    REAL(C_FLOAT), OPTIONAL, DIMENSION(*) :: Fn
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (Variant1)
+      CALL c_this__pp%Fres( T, V, n, Fres, Ft, Fv, Fn)
+    END SELECT
+  END SUBROUTINE variants__variant1_fres_wrap
+
+  SUBROUTINE variants__variant2_ctor_sub_wrap(instance, ident, nc, Tc, Vc, v&
+    &ar1, var2) BIND(C)
+    TYPE(C_PTR), VALUE :: instance
+    INTEGER(C_INT) :: ident
+    INTEGER(C_INT) :: nc
+    REAL(C_FLOAT) :: Tc
+    REAL(C_FLOAT) :: Vc
+    INTEGER(C_INT) :: var1
+    REAL(C_FLOAT) :: var2
+    TYPE(BaseEos_container_), POINTER :: instance__p
+    CALL C_F_POINTER(instance, instance__p)
+    SELECT TYPE (instance__pp => instance__p%c)
+    CLASS IS (Variant2)
+      CALL Variant2_ctor_sub(instance__pp, ident, nc, Tc, Vc, var1, var2)
+    END SELECT
+  END SUBROUTINE variants__variant2_ctor_sub_wrap
+
+  SUBROUTINE variants__variant2_db_ctor_sub_wrap(c_this, ident) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    TYPE(C_PTR), VALUE :: ident
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CHARACTER(10), POINTER :: ident__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    CALL C_F_POINTER(ident, ident__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (Variant2)
+      CALL Variant2_db_ctor_sub(c_this__pp, ident=ident__p)
+    END SELECT
+  END SUBROUTINE variants__variant2_db_ctor_sub_wrap
+
+  SUBROUTINE variants__variant2_internal_comp_wrap(c_this, T, V, n, computed&
+    &) BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: T
+    REAL(C_FLOAT) :: V
+    REAL(C_FLOAT), DIMENSION(*) :: n
+    REAL(C_FLOAT) :: computed
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (Variant2)
+      CALL c_this__pp%internal_comp( T, V, n, computed)
+    END SELECT
+  END SUBROUTINE variants__variant2_internal_comp_wrap
+
+  SUBROUTINE variants__variant2_fres_wrap(c_this, T, V, n, Fres, Ft, Fv, Fn)&
+    & BIND(C)
+    TYPE(C_PTR), VALUE :: c_this
+    REAL(C_FLOAT) :: T
+    REAL(C_FLOAT) :: V
+    REAL(C_FLOAT), DIMENSION(*) :: n
+    REAL(C_FLOAT) :: Fres
+    REAL(C_FLOAT), OPTIONAL :: Ft
+    REAL(C_FLOAT), OPTIONAL :: Fv
+    REAL(C_FLOAT), OPTIONAL, DIMENSION(*) :: Fn
+    TYPE(BaseEos_container_), POINTER :: c_this__p
+    CALL C_F_POINTER(c_this, c_this__p)
+    SELECT TYPE (c_this__pp => c_this__p%c)
+    CLASS IS (Variant2)
+      CALL c_this__pp%Fres( T, V, n, Fres, Ft, Fv, Fn)
+    END SELECT
+  END SUBROUTINE variants__variant2_fres_wrap
+
+END MODULE variants_wrap_
+
+
diff --git a/prototype/cpp_wrapper/InterfaceDefs.h b/prototype/cpp_wrapper/InterfaceDefs.h
new file mode 100644
index 00000000..1262b431
--- /dev/null
+++ b/prototype/cpp_wrapper/InterfaceDefs.h
@@ -0,0 +1,13 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#ifndef INTERFACEDEFS_H_
+#define INTERFACEDEFS_H_
+
+typedef void(*generic_fpointer)(void);
+typedef void* ADDRESS;
+
+extern "C" {
+}
+
+#endif /* INTERFACEDEFS_H_ */
diff --git a/prototype/cpp_wrapper/NotIdGas.cpp b/prototype/cpp_wrapper/NotIdGas.cpp
new file mode 100644
index 00000000..b2bf223b
--- /dev/null
+++ b/prototype/cpp_wrapper/NotIdGas.cpp
@@ -0,0 +1,22 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#include <cstring> // For strcpy
+#include "NotIdGas.h"
+
+// Constructor:
+NotIdGas::NotIdGas(int ident, int nc, float Tc, float Vc) {
+  data_ptr = nullptr;
+  allocate_notidgas_(&data_ptr); // Allocate Fortran derived type
+  notidgas_mod__notidgas_ctor_sub_wrap(data_ptr, &ident, &nc, &Tc, &Vc); // Fortran Constructor
+}
+
+// Destructor:
+NotIdGas::~NotIdGas() {
+  deallocate_notidgas_(data_ptr); // Deallocate Fortran derived type
+}
+
+void NotIdGas::Fres(float T, float V, const float n[], float* Fres, float** Ft, float** Fv, float* Fn[]) {
+  notidgas_mod__notidgas_fres_wrap(data_ptr, &T, &V, n, Fres, Ft, Fv, Fn);
+}
+
diff --git a/prototype/cpp_wrapper/NotIdGas.h b/prototype/cpp_wrapper/NotIdGas.h
new file mode 100644
index 00000000..392e924c
--- /dev/null
+++ b/prototype/cpp_wrapper/NotIdGas.h
@@ -0,0 +1,43 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#ifndef NOTIDGAS_H_
+#define NOTIDGAS_H_
+
+#ifdef _MSC_VER
+#define DLLEXPORT __declspec(dllexport)
+#else
+#define DLLEXPORT
+#endif
+
+
+#ifndef SWIG // Protect declarations from SWIG
+#include <cstdlib>
+#include "InterfaceDefs.h"
+#include "BaseEos.h"
+
+extern "C" {
+  void allocate_notidgas_(ADDRESS *caddr);
+  void deallocate_notidgas_(ADDRESS caddr);
+  void notidgas_mod__notidgas_ctor_sub_wrap(ADDRESS instance, int* ident, int* nc, float* Tc, float* Vc);
+  void notidgas_mod__notidgas_fres_wrap(ADDRESS c_this, float* T, float* V, const float n[], float* Fres, float** Ft, float** Fv, float* Fn[]);
+}
+#endif // SWIG
+
+/**
+ * \brief Not Ideal gas EoS
+ */
+class DLLEXPORT NotIdGas : public BaseEos {
+
+public:
+  NotIdGas(int ident, int nc, float Tc, float Vc);
+  ~NotIdGas();
+
+  /**
+   * \brief Compute residual and derivatives for ideal gas
+   */
+  void Fres(float T, float V, const float n[], float* Fres, float** Ft=nullptr, float** Fv=nullptr, float* Fn[]=nullptr);
+
+};
+
+#endif /* NOTIDGAS_H_ */
diff --git a/prototype/cpp_wrapper/Variant1.cpp b/prototype/cpp_wrapper/Variant1.cpp
new file mode 100644
index 00000000..87b247c6
--- /dev/null
+++ b/prototype/cpp_wrapper/Variant1.cpp
@@ -0,0 +1,39 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#include <cstring> // For strcpy
+#include "Variant1.h"
+
+// Constructor:
+Variant1::Variant1(int ident, int nc, float Tc, float Vc, int var1, float var2) {
+  data_ptr = nullptr;
+  allocate_variant1_(&data_ptr); // Allocate Fortran derived type
+  variants__variant1_ctor_sub_wrap(data_ptr, &ident, &nc, &Tc, &Vc, &var1, &var2); // Fortran Constructor
+}
+
+// Constructor:
+Variant1::Variant1(const char* ident) {
+  data_ptr = nullptr;
+  allocate_variant1_(&data_ptr); // Allocate Fortran derived type
+  // Create C array for Fortran input string data
+  char ident_c__[10+1];
+  if (ident) {
+    strncpy(ident_c__, ident, 10+1); ident_c__[10] = 0; // strncpy protects in case ident is too long
+    for (size_t i=strlen(ident_c__); i<10+1; i++) ident_c__[i] = ' '; // Add whitespace for Fortran
+  }
+  variants__variant1_db_ctor_sub_wrap(data_ptr, ident ? ident_c__ : nullptr); // Fortran Constructor
+}
+
+// Destructor:
+Variant1::~Variant1() {
+  deallocate_variant1_(data_ptr); // Deallocate Fortran derived type
+}
+
+void Variant1::internal_comp(float T, float V, const float n[], float* computed) {
+  variants__variant1_internal_comp_wrap(data_ptr, &T, &V, n, computed);
+}
+
+void Variant1::Fres(float T, float V, const float n[], float* Fres, float** Ft, float** Fv, float* Fn[]) {
+  variants__variant1_fres_wrap(data_ptr, &T, &V, n, Fres, Ft, Fv, Fn);
+}
+
diff --git a/prototype/cpp_wrapper/Variant1.h b/prototype/cpp_wrapper/Variant1.h
new file mode 100644
index 00000000..ecbf4158
--- /dev/null
+++ b/prototype/cpp_wrapper/Variant1.h
@@ -0,0 +1,43 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#ifndef VARIANT1_H_
+#define VARIANT1_H_
+
+#ifdef _MSC_VER
+#define DLLEXPORT __declspec(dllexport)
+#else
+#define DLLEXPORT
+#endif
+
+
+#ifndef SWIG // Protect declarations from SWIG
+#include <cstdlib>
+#include "InterfaceDefs.h"
+#include <string>
+#include "VariantEoS.h"
+
+extern "C" {
+  void allocate_variant1_(ADDRESS *caddr);
+  void deallocate_variant1_(ADDRESS caddr);
+  void variants__variant1_ctor_sub_wrap(ADDRESS c_this, int* ident, int* nc, float* Tc, float* Vc, int* var1, float* var2);
+  void variants__variant1_db_ctor_sub_wrap(ADDRESS c_this, const char* ident);
+  void variants__variant1_internal_comp_wrap(ADDRESS c_this, float* T, float* V, const float n[], float* computed);
+  void variants__variant1_fres_wrap(ADDRESS c_this, float* T, float* V, const float n[], float* Fres, float** Ft, float** Fv, float* Fn[]);
+}
+#endif // SWIG
+
+class DLLEXPORT Variant1 : public VariantEoS {
+
+public:
+  Variant1(int ident, int nc, float Tc, float Vc, int var1, float var2);
+  Variant1(const char* ident);
+  ~Variant1();
+
+  void internal_comp(float T, float V, const float n[], float* computed);
+
+  void Fres(float T, float V, const float n[], float* Fres, float** Ft=nullptr, float** Fv=nullptr, float* Fn[]=nullptr);
+
+};
+
+#endif /* VARIANT1_H_ */
diff --git a/prototype/cpp_wrapper/Variant2.cpp b/prototype/cpp_wrapper/Variant2.cpp
new file mode 100644
index 00000000..963f880d
--- /dev/null
+++ b/prototype/cpp_wrapper/Variant2.cpp
@@ -0,0 +1,39 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#include <cstring> // For strcpy
+#include "Variant2.h"
+
+// Constructor:
+Variant2::Variant2(int ident, int nc, float Tc, float Vc, int var1, float var2) {
+  data_ptr = nullptr;
+  allocate_variant2_(&data_ptr); // Allocate Fortran derived type
+  variants__variant2_ctor_sub_wrap(data_ptr, &ident, &nc, &Tc, &Vc, &var1, &var2); // Fortran Constructor
+}
+
+// Constructor:
+Variant2::Variant2(const char* ident) {
+  data_ptr = nullptr;
+  allocate_variant2_(&data_ptr); // Allocate Fortran derived type
+  // Create C array for Fortran input string data
+  char ident_c__[10+1];
+  if (ident) {
+    strncpy(ident_c__, ident, 10+1); ident_c__[10] = 0; // strncpy protects in case ident is too long
+    for (size_t i=strlen(ident_c__); i<10+1; i++) ident_c__[i] = ' '; // Add whitespace for Fortran
+  }
+  variants__variant2_db_ctor_sub_wrap(data_ptr, ident ? ident_c__ : nullptr); // Fortran Constructor
+}
+
+// Destructor:
+Variant2::~Variant2() {
+  deallocate_variant2_(data_ptr); // Deallocate Fortran derived type
+}
+
+void Variant2::internal_comp(float T, float V, const float n[], float* computed) {
+  variants__variant2_internal_comp_wrap(data_ptr, &T, &V, n, computed);
+}
+
+void Variant2::Fres(float T, float V, const float n[], float* Fres, float** Ft, float** Fv, float* Fn[]) {
+  variants__variant2_fres_wrap(data_ptr, &T, &V, n, Fres, Ft, Fv, Fn);
+}
+
diff --git a/prototype/cpp_wrapper/Variant2.h b/prototype/cpp_wrapper/Variant2.h
new file mode 100644
index 00000000..246e1072
--- /dev/null
+++ b/prototype/cpp_wrapper/Variant2.h
@@ -0,0 +1,43 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#ifndef VARIANT2_H_
+#define VARIANT2_H_
+
+#ifdef _MSC_VER
+#define DLLEXPORT __declspec(dllexport)
+#else
+#define DLLEXPORT
+#endif
+
+
+#ifndef SWIG // Protect declarations from SWIG
+#include <cstdlib>
+#include "InterfaceDefs.h"
+#include <string>
+#include "VariantEoS.h"
+
+extern "C" {
+  void allocate_variant2_(ADDRESS *caddr);
+  void deallocate_variant2_(ADDRESS caddr);
+  void variants__variant2_ctor_sub_wrap(ADDRESS instance, int* ident, int* nc, float* Tc, float* Vc, int* var1, float* var2);
+  void variants__variant2_db_ctor_sub_wrap(ADDRESS c_this, const char* ident);
+  void variants__variant2_internal_comp_wrap(ADDRESS c_this, float* T, float* V, const float n[], float* computed);
+  void variants__variant2_fres_wrap(ADDRESS c_this, float* T, float* V, const float n[], float* Fres, float** Ft, float** Fv, float* Fn[]);
+}
+#endif // SWIG
+
+class DLLEXPORT Variant2 : public VariantEoS {
+
+public:
+  Variant2(int ident, int nc, float Tc, float Vc, int var1, float var2);
+  Variant2(const char* ident);
+  ~Variant2();
+
+  void internal_comp(float T, float V, const float n[], float* computed);
+
+  void Fres(float T, float V, const float n[], float* Fres, float** Ft=nullptr, float** Fv=nullptr, float* Fn[]=nullptr);
+
+};
+
+#endif /* VARIANT2_H_ */
diff --git a/prototype/cpp_wrapper/VariantEoS.cpp b/prototype/cpp_wrapper/VariantEoS.cpp
new file mode 100644
index 00000000..91dd0477
--- /dev/null
+++ b/prototype/cpp_wrapper/VariantEoS.cpp
@@ -0,0 +1,10 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#include <cstring> // For strcpy
+#include "VariantEoS.h"
+
+void VariantEoS::variant_common_comp(float* mixed_common) {
+  variants__variant_common_comp_wrap(data_ptr, mixed_common);
+}
+
diff --git a/prototype/cpp_wrapper/VariantEoS.h b/prototype/cpp_wrapper/VariantEoS.h
new file mode 100644
index 00000000..8414a4a3
--- /dev/null
+++ b/prototype/cpp_wrapper/VariantEoS.h
@@ -0,0 +1,39 @@
+/* This source file automatically generated on 2023-06-07 using 
+   FortWrap wrapper generator version 2.2.2 */
+
+#ifndef VARIANTEOS_H_
+#define VARIANTEOS_H_
+
+#ifdef _MSC_VER
+#define DLLEXPORT __declspec(dllexport)
+#else
+#define DLLEXPORT
+#endif
+
+
+#ifndef SWIG // Protect declarations from SWIG
+#include <cstdlib>
+#include "InterfaceDefs.h"
+#include "BaseEos.h"
+
+extern "C" {
+  void variants__variant_common_comp_wrap(ADDRESS c_this, float* mixed_common);
+}
+#endif // SWIG
+
+class DLLEXPORT VariantEoS : public BaseEos {
+
+protected:
+  // VariantEoS can not be instantiated
+  VariantEoS() {}
+
+public:
+  virtual ~VariantEoS() {}
+
+  void variant_common_comp(float* mixed_common);
+
+  virtual void internal_comp(float T, float V, const float n[], float* computed) = 0;
+
+};
+
+#endif /* VARIANTEOS_H_ */
diff --git a/prototype/cpp_wrapper/pybind11_bindings.cpp b/prototype/cpp_wrapper/pybind11_bindings.cpp
new file mode 100644
index 00000000..6ba3beb2
--- /dev/null
+++ b/prototype/cpp_wrapper/pybind11_bindings.cpp
@@ -0,0 +1,93 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/numpy.h>
+#include <optional>
+#include "Variant2.h"
+#include "NotIdGas.h"
+#include "Variant1.h"
+namespace py = pybind11;
+using array_f = pybind11::array_t<float, pybind11::array::c_style | pybind11::array::forcecast>;
+using p_array_f = std::optional<pybind11::array_t<float, pybind11::array::c_style | pybind11::array::forcecast>>;
+
+
+#define BaseEos_bindings(Model) \
+    .def("BaseEos_init", &Model::BaseEos_init) \
+    .def("Fideal", [](Model& self, float T, float V, const std::vector<float>& n, array_f Fid, p_array_f Ft, p_array_f Fv, p_array_f Fn){ \
+            auto Fid_bind = Fid.mutable_unchecked().mutable_data(0); \
+            auto Ft_bind = (Ft.has_value()) ? Ft.value().mutable_unchecked().mutable_data(0) : nullptr; \
+            auto Fv_bind = (Fv.has_value()) ? Fv.value().mutable_unchecked().mutable_data(0) : nullptr; \
+            auto Fn_bind = (Fn.has_value()) ? Fn.value().mutable_unchecked().mutable_data(0) : nullptr; \
+            self.Fideal(T, V, n.data(), Fid_bind, &Ft_bind, &Fv_bind, &Fn_bind); \
+            }, py::arg("T")\
+            , py::arg("V")\
+            , py::arg("n")\
+            , py::arg("Fid").noconvert()\
+            , py::arg("Ft").noconvert() = py::none()\
+            , py::arg("Fv").noconvert() = py::none()\
+            , py::arg("Fn").noconvert() = py::none()\
+             \
+        ) \
+    .def("pressure", [](Model& self, float T, float V, const std::vector<float>& n, array_f p){ \
+            auto p_bind = p.mutable_unchecked().mutable_data(0); \
+            self.pressure(T, V, n.data(), p_bind); \
+            }, py::arg("T")\
+            , py::arg("V")\
+            , py::arg("n")\
+            , py::arg("p").noconvert()\
+             \
+        ) \
+    .def("set_Tc", &Model::set_Tc) \
+    .def("Fres", [](Model& self, float T, float V, const std::vector<float>& n, array_f Fres, p_array_f Ft, p_array_f Fv, p_array_f Fn){ \
+            auto Fres_bind = Fres.mutable_unchecked().mutable_data(0); \
+            auto Ft_bind = (Ft.has_value()) ? Ft.value().mutable_unchecked().mutable_data(0) : nullptr; \
+            auto Fv_bind = (Fv.has_value()) ? Fv.value().mutable_unchecked().mutable_data(0) : nullptr; \
+            auto Fn_bind = (Fn.has_value()) ? Fn.value().mutable_unchecked().mutable_data(0) : nullptr; \
+            self.Fres(T, V, n.data(), Fres_bind, &Ft_bind, &Fv_bind, &Fn_bind); \
+            }, py::arg("T")\
+            , py::arg("V")\
+            , py::arg("n")\
+            , py::arg("Fres").noconvert()\
+            , py::arg("Ft").noconvert() = py::none()\
+            , py::arg("Fv").noconvert() = py::none()\
+            , py::arg("Fn").noconvert() = py::none()\
+             \
+        ) 
+
+#define VariantEoS_bindings(Model) \
+    .def("variant_common_comp", [](Model& self, array_f mixed_common){ \
+            auto mixed_common_bind = mixed_common.mutable_unchecked().mutable_data(0); \
+            self.variant_common_comp(mixed_common_bind); \
+            }, py::arg("mixed_common").noconvert()\
+             \
+        ) \
+    .def("internal_comp", [](Model& self, float T, float V, const std::vector<float>& n, array_f computed){ \
+            auto computed_bind = computed.mutable_unchecked().mutable_data(0); \
+            self.internal_comp(T, V, n.data(), computed_bind); \
+            }, py::arg("T")\
+            , py::arg("V")\
+            , py::arg("n")\
+            , py::arg("computed").noconvert()\
+             \
+        ) 
+
+
+PYBIND11_MODULE(libdemo, handle){
+    handle.doc() = "Is this documentation? I have been told it is the best.";
+
+    py::class_<Variant2>(handle, "Variant2") 
+        .def(py::init<int, int, float, float, int, float>(), py::arg("ident"), py::arg("nc"), py::arg("Tc"), py::arg("Vc"), py::arg("var1"), py::arg("var2")) 
+        .def(py::init<const char*>(), py::arg("ident")) 
+        VariantEoS_bindings(Variant2) 
+        BaseEos_bindings(Variant2); 
+
+    py::class_<NotIdGas>(handle, "NotIdGas") 
+        .def(py::init<int, int, float, float>(), py::arg("ident"), py::arg("nc"), py::arg("Tc"), py::arg("Vc")) 
+        BaseEos_bindings(NotIdGas); 
+
+    py::class_<Variant1>(handle, "Variant1") 
+        .def(py::init<int, int, float, float, int, float>(), py::arg("ident"), py::arg("nc"), py::arg("Tc"), py::arg("Vc"), py::arg("var1"), py::arg("var2")) 
+        .def(py::init<const char*>(), py::arg("ident")) 
+        VariantEoS_bindings(Variant1) 
+        BaseEos_bindings(Variant1); 
+
+}
diff --git a/prototype/fluids/demo_db_handler.py b/prototype/fluids/demo_db_handler.py
new file mode 100644
index 00000000..f38bab40
--- /dev/null
+++ b/prototype/fluids/demo_db_handler.py
@@ -0,0 +1,179 @@
+import json
+import os
+from warnings import warn
+
+HERE = os.path.dirname(__file__)
+
+
+def find_fluid_file(ident, include_alias=False):
+    """
+    Get the absolute path to the fluid file of the fluid with identifier 'ident'
+
+    Args:
+        ident (str) : Fluid identifier
+        include_alias (bool, optional) : Also search in aliases? (Default : False)
+                                        NB: Aliases are not guaranteed to be unique!
+
+    Raises:
+        FileNotFoundError : If the fluid identifier is not found
+
+    Returns:
+        str : The absolute path to the matching fluid file
+    """
+    try:
+        filepath = f'{HERE}/{ident}.json'
+        data = json.load(open(filepath, 'r'))
+        if ident == data['ident']:
+            return filepath
+
+        if include_alias is True:
+            if ident in data['aliases']:
+                return filepath
+    except FileNotFoundError:
+        pass
+
+    for filename in os.listdir(HERE):
+        filepath = f'{HERE}/{filename}'
+        data = json.load(open(filepath, 'r'))
+        if ident == data['ident']:
+            return filepath
+
+        if include_alias is True:
+            if ident in data['aliases']:
+                return filepath
+
+    raise FileNotFoundError(f'Did not find a fluid file with identifier : {ident}')
+
+
+def add_fluid(filename, base_data):
+    """
+    Add a new fluid to the database, with minimal required data
+
+    Args:
+        filename (str) : Name of the fluid file (including file ending, which should be .json)
+        base_data (dict) : A dict matching 'default_base_data', which supplies values for all parameters that are 'None'
+                            in 'default_base_data'.
+    """
+    filepath = f'{HERE}/{filename}'
+    if os.path.isfile(filepath):
+        raise FileExistsError(f'A fluid file already exists at {filepath}')
+
+    base_data = verify_and_amend_dataset('base_model', base_data)
+
+    if base_data['ident'] != filename.split('.')[0]:
+        warn(f'Fluid identifier {base_data["ident"]} does not match filename {filename}', SyntaxWarning, stacklevel=2)
+
+    json.dump(base_data, open(filepath, 'w'), indent=4)
+    print(f'Created Fluid file {filepath}')
+
+
+def verify_and_amend_dataset(model, model_data):
+    """
+    Check that the dict 'model_data' contains all required parameters. If optional parameters are omitted, insert
+    the default values.
+
+    Args:
+        model (str) : The model key
+        model_data (dict) : The model parameter set
+
+    Raises:
+        KeyError : If 'model' does not have a default data set to verify against
+        KeyError : If 'model_data' is missing required parameters
+        SyntaxWarning : If 'model_data' contains parameters (other than "note", "comment" or "doi") that are not found
+                        in the default parameter set.
+
+    Returns:
+        dict : The supplied 'model_data' with inserted default values if optional values were not supplied
+    """
+
+    if model == 'base_model':
+        default_data = {'ident': None,  # Fluid identifier (str)
+                        'dipole_moment' : 0, # float, optional
+                        'T_crit': None,  # float
+                        'T_triple': None,  # float
+                        'p_triple': None,  # float
+                        }
+    elif model in ('Variant1', 'Variant2'):
+        default_data = {'p1' : None,
+                        'p2' : None}
+    elif model == 'NotIdGas':
+        default_data = {}
+    elif model == 'CritModel':
+        default_data = {'T_crit' : None,
+                        'p_crit' : None}
+    else:
+        raise KeyError(f'Unknown model type {model}')
+
+    for k in model_data.keys():
+        if (k not in default_data.keys()) and (k.lower() not in ('note', 'comment', 'doi', 'alt_reference')):
+            warn(f'A value for {k} was supplied when creating Model type : {model}, '
+                 f'but no corresponding entry exists in the default set!', SyntaxWarning, stacklevel=2)
+
+    for k, v in default_data.items():
+        if (v is None) and k not in model_data.keys():
+            raise KeyError(f'A value for {k} must be supplied to create {model}')
+        elif k not in model_data.keys():
+            model_data[k] = v
+
+    return model_data
+
+
+def add_model_to_fluid(ident, model, reference, model_data, make_default=False):
+    """
+    Add a new model parameter set for the fluid 'ident'
+
+    Args:
+        ident (str) : fluid identifier (not alias)
+        model (str) : Model identifier
+        reference (str) : Parameter reference
+        model_data (dict) : Model parameters
+        make_default (bool, optional) : Make this parameter set the default for this model? (Default : False)
+    """
+    filepath = find_fluid_file(ident)
+    data = json.load(open(filepath, 'r'))
+
+    if model not in data.keys():
+        data[model] = {}
+        make_default = True
+
+    if 'alt_reference' in model_data.keys():
+        alt_reference = model_data['alt_reference']
+    else:
+        alt_reference = None
+
+    if reference.lower() != 'default':
+        if reference in data[model].keys():
+            raise KeyError(f'Parameter reference : {reference} already exists for Fluid {ident}, model {model}')
+        if 'default' in data[model].keys():
+            if reference == data[model]['default']['alt_reference']:
+                raise KeyError(f'Parameter reference : {reference} already exists for Fluid {ident}, model {model}')
+        if make_default is True:
+            model_data['alt_reference'] = reference
+            add_model_to_fluid(ident, model, 'default', model_data)
+            return
+    else:
+        if alt_reference is None:
+            raise KeyError('An alternative reference must be supplied when setting the default parameter set.')
+        if alt_reference in data[model].keys():
+            raise KeyError(f'The alternative reference {alt_reference} already exists for Fluid {ident}, model {model}')
+        if 'default' in data[model].keys():
+            if alt_reference == data[model]['default']['alt_reference']:
+                raise KeyError(f'The alternative reference {alt_reference} already exists for '
+                               f'Fluid {ident}, model {model}')
+
+    if (reference.lower() == 'default') and ('default' in data[model].keys()):
+        old_default = data[model]['default'].copy()
+        old_default.pop('alt_reference')
+        data[model][data[model]['default']['alt_reference']] = old_default
+
+    model_data = verify_and_amend_dataset(model, model_data)
+    data[model][reference] = model_data
+    json.dump(data, open(filepath, 'w'), indent=4)
+    print(f'Added parameter set {reference : <20} {model : <15} {ident : <10}'
+          + (f' (Alternative ref : {alt_reference})' if alt_reference is not None else ''))
+
+
+if __name__ == '__main__':
+    sA = {'T_crit' : 105, 'p_crit' : 15}
+    sA_2 = {'T_crit': 110, 'p_crit' : 20}
+    add_model_to_fluid('speciesA', 'CritModel', 'Betterref', sA_2, make_default=True)
\ No newline at end of file
diff --git a/prototype/fluids/fluid_db_handler.py b/prototype/fluids/fluid_db_handler.py
new file mode 100644
index 00000000..5ccfa328
--- /dev/null
+++ b/prototype/fluids/fluid_db_handler.py
@@ -0,0 +1,177 @@
+import json
+import os
+from warnings import warn
+
+HERE = os.path.dirname(__file__)
+
+def find_fluid_file(ident, include_alias=False):
+    """
+    Get the absolute path to the fluid file of the fluid with identifier 'ident'
+
+    Args:
+        ident (str) : Fluid identifier
+        include_alias (bool, optional) : Also search in aliases? (Default : False)
+                                        NB: Aliases are not guaranteed to be unique!
+
+    Raises:
+        FileNotFoundError : If the fluid identifier is not found
+
+    Returns:
+        str : The absolute path to the matching fluid file
+    """
+    try:
+        filepath = f'{HERE}/{ident}.json'
+        data = json.load(open(filepath, 'r'))
+        if ident == data['ident']:
+            return filepath
+
+        if include_alias is True:
+            if ident in data['aliases']:
+                return filepath
+    except FileNotFoundError:
+        pass
+
+    for filename in os.listdir(HERE):
+        filepath = f'{HERE}/{filename}'
+        data = json.load(open(filepath, 'r'))
+        if ident == data['ident']:
+            return filepath
+
+        if include_alias is True:
+            if ident in data['aliases']:
+                return filepath
+
+    raise FileNotFoundError(f'Did not find a fluid file with identifier : {ident}')
+
+def add_fluid(filename, base_data):
+    """
+    Add a new fluid to the database, with minimal required data
+
+    Args:
+        filename (str) : Name of the fluid file (including file ending, which should be .json)
+        base_data (dict) : A dict matching 'default_base_data', which supplies values for all parameters that are 'None'
+                            in 'default_base_data'.
+    """
+    filepath = f'{HERE}/{filename}'
+    if os.path.isfile(filepath):
+        raise FileExistsError(f'A fluid file already exists at {filepath}')
+
+    base_data = verify_and_amend_dataset('base_model', base_data)
+
+    if base_data['ident'] != filename.split('.')[0]:
+        warn(f'Fluid identifier {base_data["ident"]} does not match filename {filename}', SyntaxWarning, stacklevel=2)
+
+    json.dump(base_data, open(filepath, 'w'), indent=4)
+    print(f'Created Fluid file {filepath}')
+
+def verify_and_amend_dataset(model, model_data):
+    """
+    Check that the dict 'model_data' contains all required parameters. If optional parameters are omitted, insert
+    the default values.
+
+    Args:
+        model (str) : The model key
+        model_data (dict) : The model parameter set
+
+    Raises:
+        KeyError : If 'model' does not have a default data set to verify against
+        KeyError : If 'model_data' is missing required parameters
+        SyntaxWarning : If 'model_data' contains parameters (other than "note", "comment" or "doi") that are not found
+                        in the default parameter set.
+
+    Returns:
+        dict : The supplied 'model_data' with inserted default values if optional values were not supplied
+    """
+
+    if model == 'base_model':
+        default_data = {'ident': None,  # Fluid identifier (str)
+                         'aliases': (),  # Alternative fluid identifiers (tuple[str]) (optional)
+                         'CAS_number': None,  # int
+                         'SMILES': None,  # str
+                         'mole_weight': None,  # float
+                         'T_crit': None,  # float
+                         'p_crit': None,  # float
+                         'rho_crit': None,  # float
+                         'acentric_factor': None,  # float
+                         'T_triple': None,  # float
+                         'p_triple': None,  # float
+                         'rho_triple': None,  # float
+                         'melting_enthalpy_triple': None,  # float
+                         'melting_entropy_triple': None,  # float
+                         'dipole_moment': 0,  # float (optional)
+                         'quadrupole_moment': 0,  # float (optional
+                         'formation_enthalpy': None,  # float
+                         'formation_entropy': None,  # float
+                         'default_ideal_cp_correlation': None,  # int
+                         'default_ideal_cp_coefficients': None  # tuple[float]
+                         }
+    elif model == 'SAFTVRMIE':
+        default_data = {'m' : None,
+                        'eps_div_k' : None,
+                        'sigma' : None,
+                        'lambda_a' : None,
+                        'lambda_r' : None,
+                        }
+    else:
+        raise KeyError(f'Unknown model type {model}')
+
+    for k in model_data.keys():
+        if (k not in default_data.keys()) and (k.lower() not in ('note', 'comment', 'doi')):
+            warn(f'A value for {k} was supplied when creating Model type : {model}, '
+                 f'but no corresponding entry exists in the default set!', SyntaxWarning, stacklevel=2)
+
+    for k, v in default_data.items():
+        if (v is None) and k not in model_data.keys():
+            raise KeyError(f'A value for {k} must be supplied to create {model}')
+        elif k not in model_data.keys:
+            model_data[k] = v
+
+    return model_data
+
+def add_model_to_fluid(ident, model, reference, model_data, make_default=False):
+    """
+    Add a new model parameter set for the fluid 'ident'
+
+    Args:
+        ident (str) : fluid identifier (not alias)
+        model (str) : Model identifier
+        reference (str) : Parameter reference
+        model_data (dict) : Model parameters
+        make_default (bool, optional) : Make this parameter set the default for this model? (Default : False)
+    """
+    filepath = find_fluid_file(ident)
+    data = json.load(open(filepath, 'r'))
+
+    if model not in data.keys():
+        data[model] = {}
+        make_default = True
+
+    if 'alt_reference' in model_data.keys():
+        alt_reference = model_data['alt_reference']
+    else:
+        alt_reference = None
+
+    if reference.lower() != 'default':
+        if (reference in data[model].keys()) or (reference == data[model]['default']['alt_reference']):
+            raise KeyError(f'Parameter reference : {reference} already exists for Fluid {ident}, model {model}')
+        if make_default is True:
+            model_data['alt_reference'] = reference
+            add_model_to_fluid(ident, model, 'default', model_data)
+    else:
+        if alt_reference is None:
+            raise KeyError('An alternative reference must be supplied when setting the default parameter set.')
+        elif (alt_reference in data[model].keys()) or (alt_reference == data[model]['default']['alt_reference']):
+            raise KeyError(f'The alternative reference {alt_reference} already exists for Fluid {ident}, model {model}')
+
+
+    if (reference.lower() == 'default') and ('default' in data[model].keys()):
+        data[model][data[model]['defalut']['alt_reference']] = data[model]['default']
+
+    model_data = verify_and_amend_dataset(model, model_data)
+    data[model][reference] = model_data
+    json.dump(data, open(filepath, 'w'), indent=4)
+    print(f'Added parameter set {reference} for {model} to {ident}')
+
+if __name__ == '__main__':
+
+    add_fluid('speciesA', )
\ No newline at end of file
diff --git a/prototype/fluids/speciesA.json b/prototype/fluids/speciesA.json
new file mode 100644
index 00000000..f506c6a8
--- /dev/null
+++ b/prototype/fluids/speciesA.json
@@ -0,0 +1,48 @@
+{
+    "ident": "speciesA",
+    "T_crit": 100,
+    "T_triple": 20,
+    "p_triple": 150,
+    "dipole_moment": 0,
+    "Variant1": {
+        "default": {
+            "p1": 30,
+            "p2": 150,
+            "alt_reference": "Best"
+        },
+        "otherdummy": {
+            "p1": 50,
+            "p2": 15
+        },
+        "Dummyref": {
+            "p1": 10,
+            "p2": 20
+        }
+    },
+    "Variant2": {
+        "default": {
+            "p1": -30,
+            "p2": -150,
+            "alt_reference": "Best"
+        },
+        "otherdummy": {
+            "p1": -50,
+            "p2": -15
+        },
+        "Dummyref": {
+            "p1": -10,
+            "p2": -20
+        }
+    },
+    "CritModel": {
+        "default": {
+            "T_crit": 110,
+            "p_crit": 20,
+            "alt_reference": "Betterref"
+        },
+        "Dummyref": {
+            "T_crit": 105,
+            "p_crit": 15
+        }
+    }
+}
\ No newline at end of file
diff --git a/prototype/fluids/speciesB.json b/prototype/fluids/speciesB.json
new file mode 100644
index 00000000..6c8f7ab3
--- /dev/null
+++ b/prototype/fluids/speciesB.json
@@ -0,0 +1,37 @@
+{
+    "ident": "speciesB",
+    "dipole_moment": 10,
+    "T_crit": 500,
+    "T_triple": 200,
+    "p_triple": 15,
+    "Variant1": {
+        "default": {
+            "p1": 45,
+            "p2": 165,
+            "alt_reference": "Best"
+        },
+        "otherdummy": {
+            "p1": 65,
+            "p2": 30
+        },
+        "Dummyref": {
+            "p1": 25,
+            "p2": 35
+        }
+    },
+    "Variant2": {
+        "default": {
+            "p1": -20,
+            "p2": -140,
+            "alt_reference": "Best"
+        },
+        "otherdummy": {
+            "p1": -40,
+            "p2": -5
+        },
+        "Dummyref": {
+            "p1": 0,
+            "p2": -10
+        }
+    }
+}
\ No newline at end of file
diff --git a/prototype/py_wrapper/BaseEoS.py b/prototype/py_wrapper/BaseEoS.py
new file mode 100644
index 00000000..be9f02c6
--- /dev/null
+++ b/prototype/py_wrapper/BaseEoS.py
@@ -0,0 +1,54 @@
+from . import libdemo # The module defined by pybind11
+import numpy as np
+from collections.abc import Iterable
+
+class BaseEoS:
+
+    def __init__(self):
+        self.nc = 0
+        self.lib = None # In this variable, we will store an initialized C++-side model.
+
+    def pressure(self, T, V, n):
+        p = np.array([0], dtype=np.float32)
+        self.lib.pressure(T, V, n, p)
+        return p[0]
+
+    def get_empty_arrays(self, flags, dims, dtypes):
+        """
+        Helper function to generate empty arrays that are used to pass optional, modifiable arguments to the C++ side.
+
+        Note that optional floats that are to be modified must be passed as 1d arrays with a single element.
+        """
+        if not isinstance(dims, Iterable):
+            dims = dims * np.ones_like(flags)
+        if not isinstance(dtypes, Iterable):
+            dtypes = [dtypes for _ in flags]
+
+        arrs = [np.empty(tuple(dims[fi] for _ in range(dims[fi])), dtype=dtypes[fi]) if (flags[fi] is not None) else
+                None
+                for fi in range(len(flags))]
+        return arrs
+
+    def Fres(self, T, V, n, Ft=None, Fv=None, Fn=None):
+        F = np.empty((1,), dtype=np.float32)
+        flags = [Ft, Fv, Fn]
+        dims = [1, 1, 2]
+        dtypes = np.float32
+        # These arrays will be modified. If they contain None, Pybind will ensure that nullptr is passed to C++,
+        # Such that Fortran never sees them.
+        Ft, Fv, Fn = self.get_empty_arrays(flags, dims, dtypes)
+
+        self.lib.Fres(T, V, n, F, Ft, Fv, Fn)
+
+        return F[0], [i for i in Ft] if (Ft is not None) else None, \
+                    [i for i in Fv] if (Fv is not None) else None, \
+                    [i for i in Fn] if (Fn is not None) else None
+
+class NotIdGas(BaseEoS):
+
+    def __init__(self, ident, ncomps, Tc, Vc):
+        self.nc = ncomps
+        super().__init__()
+        self.lib = libdemo.NotIdGas(ident, ncomps, Tc, Vc) # Initialize the corresponding C++ side model
+
+
diff --git a/prototype/rebuild_and_run.sh b/prototype/rebuild_and_run.sh
new file mode 100644
index 00000000..d1bd6010
--- /dev/null
+++ b/prototype/rebuild_and_run.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+set -e
+
+g++-12 -c main.cpp
+g++-12 -o main main.o wrappers/*.o src/*.o /opt/homebrew/lib/gcc/12/libgfortran.dylib
+./main
\ No newline at end of file
diff --git a/prototype/rebuild_fortran_prog.sh b/prototype/rebuild_fortran_prog.sh
new file mode 100644
index 00000000..6f31f826
--- /dev/null
+++ b/prototype/rebuild_fortran_prog.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+set -e
+
+cd src
+
+rm ./*.o
+rm ./*.mod
+
+gfortran -c compdata.f90 compdatadb.f90 compdatautils.f90 constants.f90 base_eos.f90 variants.f90 notidgas.f90
+
+cp constants.mod ../fortran_test/constants.mod
+cp base_eos.mod ../fortran_test/base_eos.mod
+cp notidgas_mod.mod ../fortran_test/notidgas_mod.mod
+cp variants.mod ../fortran_test/variants.mod
+cp compdata_mod.mod ../fortran_test/compdata_mod.mod
+cp compdatadb.mod ../fortran_test/compdatadb.mod
+cp compdatautils.mod ../fortran_test/compdatautils.mod
+
+cd ../fortran_test
+
+gfortran -c main.f90
+gfortran -o fortran_main main.o ../src/*.o
+./fortran_main
\ No newline at end of file
diff --git a/prototype/src/CMakeLists.txt b/prototype/src/CMakeLists.txt
new file mode 100644
index 00000000..f8eadb43
--- /dev/null
+++ b/prototype/src/CMakeLists.txt
@@ -0,0 +1,2 @@
+file(GLOB libsources *.f90)
+add_library(demopack ${libsources})
\ No newline at end of file
diff --git a/prototype/src/base_eos.f90 b/prototype/src/base_eos.f90
new file mode 100644
index 00000000..93a6fdb5
--- /dev/null
+++ b/prototype/src/base_eos.f90
@@ -0,0 +1,98 @@
+MODULE base_eos
+  use constants, only: Rgas
+  use compdata_mod, only: BaseData
+  use compdatautils, only: get_base_data
+
+  IMPLICIT NONE
+  
+  PRIVATE
+  PUBLIC :: BaseEos, BaseEoS_init
+
+ !> Base EoS class
+  TYPE, ABSTRACT :: BaseEos
+    INTEGER :: ident
+    integer :: ncomps 
+    real :: Tc, Vc, Tt, pt, dipole
+  CONTAINS
+    ! All inherriting classes must implement this
+    PROCEDURE (Fres_template), DEFERRED :: Fres
+    ! These are implemented at the top level
+    PROCEDURE :: Fideal, pressure
+    !> A set function - also implemented at the top level
+    PROCEDURE :: set_Tc
+  END TYPE BaseEos
+  
+  ABSTRACT INTERFACE
+    !> Residual helmholtz energy
+    subroutine Fres_template(this, T, V, n, Fres, Ft, Fv, Fn)
+      IMPORT :: BaseEos
+      CLASS(BaseEos), INTENT(in) :: this
+      real, intent(in) :: T, V
+      real, dimension(this%ncomps), intent(in) :: n
+      real, intent(out) :: Fres
+      real, optional, intent(out) :: Ft, Fv
+      real, optional, dimension(this%ncomps), intent(out) :: Fn
+    END subroutine Fres_template
+  END INTERFACE
+
+CONTAINS
+
+subroutine BaseEos_init(this, ident)
+    class(BaseEoS), intent(inout) :: this
+    character(len=10), intent(in) :: ident
+    type(BaseData) :: bdata
+    bdata = get_base_data(ident)
+
+    this%Tc = bdata%Tc
+    this%Vc = bdata%Vc
+    this%Tt = bdata%Tt
+    this%pt = bdata%pt
+    this%dipole = bdata%dipole
+end subroutine BaseEos_init
+
+subroutine Fideal(this, T, V, n, Fid, Ft, Fv, Fn)
+    CLASS(BaseEos), INTENT(in) :: this
+    real, intent(in) :: T, V
+    real, dimension(this%ncomps), intent(in) :: n
+    real, intent(out) :: Fid
+    real, optional, intent(out) :: Ft, Fv
+    real, optional, dimension(this%ncomps), intent(out) :: Fn
+    real :: V0 = 0.25
+    real :: T0 = 300.0
+    integer :: i
+
+    Fid = sum(n) * 3.0 * Rgas * log(T / T0) / 2.0 + sum(n) * Rgas * T * log(V0 / V)
+    if (present(Ft)) Ft = sum(n) * 3.0 * Rgas / 2.0 
+    if (present(Fv)) Fv = - sum(n) * Rgas * T / V
+    if (present(Fn)) then
+        do i = 1, this%ncomps
+            Fn(i) = 3.0 * Rgas * T / 2.0 - Rgas * T / V
+        enddo
+    endif
+
+END subroutine Fideal
+
+  subroutine pressure(this, T, V, n, p)
+    CLASS(BaseEos), INTENT(in) :: this
+    real, intent(in) :: T, V
+    real, dimension(this%ncomps), intent(in) :: n
+    real, intent(out) :: p
+    real :: F, Ft
+    real, dimension(this%ncomps) :: Fn
+    real :: Fv_res, Fv_id
+
+    call this%Fideal(T, V, n, F, Ft, Fv_id, Fn)
+    call this%Fres(T, V, n, F, Ft, Fv_res, Fn)
+
+    p = - (Fv_res + Fv_id)
+
+  END subroutine pressure
+
+  subroutine set_Tc(this, Tc)
+    class(BaseEos), intent(inout) :: this
+    real, intent(in) :: Tc
+
+    this%Tc = Tc
+  end subroutine set_Tc
+
+END MODULE base_eos
diff --git a/prototype/src/compdata.f90 b/prototype/src/compdata.f90
new file mode 100644
index 00000000..ac932144
--- /dev/null
+++ b/prototype/src/compdata.f90
@@ -0,0 +1,31 @@
+module compdata_mod
+
+implicit none
+
+! Struct to hold the data all fluids require
+type BaseData
+    real :: Tc
+    real :: Vc
+    real :: Tt
+    real :: pt
+    real :: dipole
+end type BaseData
+
+! Struct to hold parameters for VariantEoS models
+type VariantData
+    character(len=10) :: ref
+    real :: p1
+    real :: p2
+end type VariantData
+
+! Top level database structure
+! The fluid database consists of an array of CompData objects, found in compdatadb
+! To retrieve the compdata object corresponding to a specific identifier, use the methods found there
+type CompData
+    character(len=10) :: ident
+
+    type(BaseData) :: base
+    type(VariantData), allocatable, dimension(:) :: VariantDatadb
+end type CompData
+
+end module compdata_mod
\ No newline at end of file
diff --git a/prototype/src/compdatadb.f90 b/prototype/src/compdatadb.f90
new file mode 100644
index 00000000..b6cc10ad
--- /dev/null
+++ b/prototype/src/compdatadb.f90
@@ -0,0 +1,91 @@
+! This is the module that holds the entire fluid database in the form of the "datadb" array
+! The array should not be accessed directly, but through the get_compdata subroutine.
+! This file is intended to be automatically generated from the json files in the 'fluids' directory.
+
+module compdatadb
+use compdata_mod
+implicit none
+
+public :: datadb_len, get_compdata, str_eq
+
+integer :: datadb_len=2
+logical, save :: datadb_is_initialized = .false.
+type(CompData), dimension(2), save :: datadb
+
+contains
+
+! Method to retrieve the CompData object corresponding to a specific fluid
+! Ensures that the compdatadb array is only initialized once.
+subroutine get_compdata(ident, data)
+    character(len=10), intent(in) :: ident
+    type(Compdata), intent(out) :: data
+    integer :: i
+
+    if (.not. datadb_is_initialized) then
+        call init_datadb()
+    endif
+
+    do i = 1, datadb_len
+        if (str_eq(ident, datadb(i)%ident)) then
+            data = datadb(i)
+            return
+        endif
+    enddo
+
+    print*, "No data found for ", ident
+    error stop
+end subroutine get_compdata
+
+function str_eq(str1, str2) result(equal)
+    character(len=*), intent(in) :: str1, str2
+    logical :: equal
+    integer :: i
+
+    equal = .true.
+
+    do i = 1, len(str1)
+        if (str1(i:i) /= str2(i:i)) then
+            equal = .false.
+            exit
+        end if
+        if (str1(i:i) == char(0)) exit
+    end do
+
+    if (len(str1) /= len(str2)) equal = .false.
+
+end function str_eq
+
+subroutine init_datadb()
+datadb(1)%ident = "speciesA" // char(0)
+    datadb(1)%base%Tc = 100
+    datadb(1)%base%Vc = 10
+    datadb(1)%base%Tt = 30
+    datadb(1)%base%pt = 15
+    allocate(datadb(1)%variantdatadb(2))
+        datadb(1)%variantdatadb(1)%ref = "default" // char(0)
+        datadb(1)%variantdatadb(1)%p1 = 10
+        datadb(1)%variantdatadb(1)%p2 = 25
+
+        datadb(1)%variantdatadb(2)%ref = "ref2" // char(0)
+        datadb(1)%variantdatadb(2)%p1 = 100
+        datadb(1)%variantdatadb(2)%p2 = 250
+
+datadb(2)%ident = "speciesB" // char(0)
+    datadb(2)%base%Tc = 500
+    datadb(2)%base%Vc = 30
+    datadb(2)%base%Tt = 300
+    datadb(2)%base%pt = 150
+
+    allocate(datadb(2)%variantdatadb(2))
+        datadb(2)%variantdatadb(1)%ref = "default" // char(0)
+        datadb(2)%variantdatadb(1)%p1 = -10
+        datadb(2)%variantdatadb(1)%p2 = -25
+
+        datadb(2)%variantdatadb(2)%ref = "ref2" // char(0)
+        datadb(2)%variantdatadb(2)%p1 = -100
+        datadb(2)%variantdatadb(2)%p2 = -250
+
+datadb_is_initialized = .true.
+end subroutine init_datadb
+
+end module compdatadb
\ No newline at end of file
diff --git a/prototype/src/compdatautils.f90 b/prototype/src/compdatautils.f90
new file mode 100644
index 00000000..ed52303f
--- /dev/null
+++ b/prototype/src/compdatautils.f90
@@ -0,0 +1,49 @@
+! Functions and subroutines for accessing the datadb array in compdatadb.f90
+! Specifically, to retrieve the BaseData struct for a specific fluid, and the VariantData struct for a
+! Specific fluid and parameter reference
+module compdatautils
+
+use compdata_mod, only : CompData, VariantData, BaseData
+use compdatadb, only : get_compdata, datadb_len, str_eq
+implicit none
+
+private
+public :: get_base_data, get_variant_compdata
+
+contains
+
+! Get the BaseData struct corresponding to a specific fluid
+function get_base_data(ident) result(bdata)
+    character(len=10), intent(in) :: ident
+    type(BaseData) :: bdata
+    type(CompData) :: cdata
+
+    call get_compdata(ident, cdata)
+    bdata = cdata%base
+end function get_base_data
+
+! Get the VariantData struct corresponding to a specific fluid and parameter reference
+function get_variant_compdata(ident, ref) result(data)
+    character(len=10), intent(in) :: ident
+    character(len=10), intent(in) :: ref
+    type(VariantData) :: data
+
+    type(CompData) :: cdata
+    integer :: i
+
+    call get_compdata(ident, cdata)
+
+    print*, "Searching for ref : ", ref
+    do i = 1, size(cdata%VariantDataDb)
+        print*, "Checking : ", cdata%VariantDataDb(i)%ref
+        if (str_eq(ref, cdata%VariantDataDb(i)%ref)) then
+            data = cdata%VariantDataDb(i)
+            return
+        endif
+    enddo
+
+    print*, "No data found for ", ident, ", ref : ", ref
+    error stop
+end function get_variant_compdata
+
+end module compdatautils
\ No newline at end of file
diff --git a/prototype/src/constants.f90 b/prototype/src/constants.f90
new file mode 100644
index 00000000..35eac1d2
--- /dev/null
+++ b/prototype/src/constants.f90
@@ -0,0 +1,10 @@
+module constants
+
+implicit NONE
+private 
+
+public Rgas
+
+real :: Rgas = 8.314
+
+end module constants
\ No newline at end of file
diff --git a/prototype/src/notidgas.f90 b/prototype/src/notidgas.f90
new file mode 100644
index 00000000..88449a9a
--- /dev/null
+++ b/prototype/src/notidgas.f90
@@ -0,0 +1,59 @@
+MODULE notidgas_mod
+  use base_eos, only: BaseEos
+  IMPLICIT NONE
+  
+  PRIVATE
+  PUBLIC :: NotIdGas, NotIdGas_ctor_sub
+
+  !> Not Ideal gas EoS
+  TYPE, EXTENDS(BaseEos) :: NotIdGas
+  CONTAINS
+    PROCEDURE :: Fres => NotIdGas_Fres
+  END TYPE NotIdGas
+
+  INTERFACE NotIdGas
+    PROCEDURE NotIdGas_ctor
+  END INTERFACE NotIdGas
+
+Contains
+
+  function NotIdGas_ctor(ident, ncomps, Tc, Vc) result(instance)
+    INTEGER, INTENT(in) :: ident
+    integer, intent(in) :: ncomps
+    real, intent(in) :: Tc, Vc ! Read these from some database
+    TYPE (NotIdGas) :: instance
+
+    instance%ident = ident
+    instance%ncomps = ncomps
+    instance%Tc = Tc
+    instance%Vc = Vc
+  END function NotIdGas_ctor
+
+  subroutine NotIdGas_ctor_sub(instance, ident, nc, Tc, Vc)
+    TYPE (NotIdGas) :: instance
+    INTEGER, INTENT(in) :: ident
+    integer, intent(in) :: nc
+    real, intent(in) :: Tc, Vc ! Read these from some database
+    
+    instance = NotIdGas(ident, nc, Tc, Vc)
+  END subroutine NotIdGas_ctor_sub
+
+  !> Compute residual and derivatives for ideal gas
+  subroutine NotIdGas_Fres(this, T, V, n, Fres, Ft, Fv, Fn)
+    CLASS(NotIdGas), INTENT(in) :: this
+    real, intent(in) :: T, V
+    real, dimension(this%ncomps), intent(in) :: n
+    real, intent(out) :: Fres
+    real, optional, intent(out) :: Ft, Fv
+    real, optional, dimension(this%ncomps), intent(out) :: Fn
+    integer :: i
+
+    Fres = T / this%Tc + (V / sum(n)) / this%Vc
+    Ft = 1 / this%Tc
+    Fv = 1 / (sum(n) * this%Vc) 
+    do i = 1,this%ncomps
+        Fn(i) = 1 / (this%Vc * n(i))
+    enddo
+  END subroutine NotIdGas_Fres
+
+end MODULE notidgas_mod
\ No newline at end of file
diff --git a/prototype/src/variants.f90 b/prototype/src/variants.f90
new file mode 100644
index 00000000..586c0c61
--- /dev/null
+++ b/prototype/src/variants.f90
@@ -0,0 +1,263 @@
+! Module demonstrating chained inheritance
+! The class VariantEoS is a dummy for e.g. ''cubic'', 'saft' etc.
+!   It contains:
+!   Two Parameters (p1, p2) not contained in BaseEoS, that each inheriting EoS must set.
+!   One procedure that is common to all inheriting classes (implemented)
+!       returns: The average of the two parameters (p1 + p2) / 2.0
+!   Procedure internal_comp (deferred)
+
+! Variant1 
+!   Implements internalcomp
+!     Returns: T * p1 + V * p2
+!   Fres:
+!     Returns:
+!       Fres: ((p1 + p2) / 2.0) + T * p1 + V * p2
+!       Fv = Ft = Fn = 0.0
+
+! Variant2
+!   Implements internalcomp
+!     Returns: V * p1 - T * p2
+!   Fres:
+!     Returns:
+!       Fres: ((p1 + p2) / 2.0) - (V * p1 - T * p2)
+!       Fv : sum(n) * Rgas * T / V
+!       Ft = Fn = 0.0
+
+module variants
+  use base_eos, only: BaseEos, BaseEos_init
+  use constants, only: Rgas
+  use compdata_mod, only: VariantData
+  use compdatautils, only: get_variant_compdata
+
+  implicit none
+  private
+  public :: VariantEoS, Variant1, Variant2, &
+            Variant1_ctor_sub, Variant1_db_ctor_sub, Variant1_db_ctor, &
+            Variant2_ctor_sub, Variant2_db_ctor_sub, Variant2_db_ctor
+
+  ! Intermediate class for many EoS's that are similar,
+  ! e.g. cubics, saft-type, etc.
+  type, abstract, extends(BaseEos) :: VariantEoS
+    ! Contains some parameters that are common to e.g. all cubics
+    integer :: param1
+    real :: param2 
+  contains
+    ! Some internal computation that is common to e.g. all cubics
+    ! For example mixing rules.
+    procedure :: variant_common_comp
+
+    ! Some method that all inherriting classes must implement
+    procedure (internal_comp_template), deferred :: internal_comp
+  end type VariantEoS
+
+  abstract interface
+    subroutine internal_comp_template(this, T, V, n, computed)
+      import VariantEoS
+      class(VariantEoS), intent(in) :: this
+      real, intent(in) :: T, V
+      real, dimension(this%ncomps), intent(in) :: n
+      real, intent(out) :: computed
+    end subroutine internal_comp_template
+  end interface
+
+  type, extends(VariantEoS) :: Variant1
+    integer :: variant1_specific_param
+  contains
+    procedure :: internal_comp => variant1_internal_comp
+    procedure :: Fres => variant1_Fres
+  end type Variant1
+
+  interface Variant1
+    procedure Variant1_ctor
+    procedure Variant1_db_ctor
+  end interface Variant1
+
+  type, extends(VariantEoS) :: Variant2
+    integer :: variant2_specific_param
+  contains
+    procedure :: internal_comp => variant2_internal_comp
+    procedure :: Fres => variant2_Fres
+  end type Variant2
+
+  interface Variant2
+    procedure Variant2_ctor
+    procedure Variant2_db_ctor
+  end interface Variant2
+
+contains
+
+subroutine variant_common_comp(this, mixed_common)
+    class(VariantEoS), intent(in) :: this
+    real, intent(out) :: mixed_common
+
+    mixed_common = (this%param1 + this%param2) / 2.0
+end subroutine variant_common_comp
+
+function Variant1_ctor(ident, nc, Tc, Vc, var1, var2) result(instance)
+  integer, intent(in) :: ident
+  integer, intent(in) :: nc
+  real, intent(in) :: Tc, Vc
+  integer, intent(in) :: var1
+  real, intent(in) :: var2
+  type(Variant1) :: instance
+
+  instance%ident = ident
+  instance%ncomps = nc
+  instance%Tc = Tc
+  instance%Vc = Vc
+  instance%param1 = var1
+  instance%param2 = var2
+end function Variant1_ctor
+
+subroutine Variant1_ctor_sub(this, ident, nc, Tc, Vc, var1, var2)
+    TYPE (Variant1) :: this
+    INTEGER, INTENT(in) :: ident
+    integer, intent(in) :: nc
+    real, intent(in) :: Tc, Vc ! Read these from some database
+    integer, intent(in) :: var1
+    real, intent(in) :: var2
+    this = Variant1(ident, nc, Tc, Vc, var1, var2)
+end subroutine Variant1_ctor_sub
+
+function Variant1_db_ctor(ident, ref) result(instance)
+    character(len=10), intent(in) :: ident
+    character(len=10), optional, intent(in) :: ref
+    character(len=10) :: ref_ = "default" // char(0)
+    type(Variant1) :: instance
+    type(VariantData) :: data
+
+    if (present(ref)) then
+        ref_ = ref
+    end if
+
+    print*, "Getting VariantData for : ", ident, ", ", ref_
+
+    data = get_variant_compdata(ident, ref_)
+
+    print*, "Got VariantData for : ", ident, ", ", ref_
+
+    call BaseEos_init(instance, ident)
+
+    instance%param1 = data%p1
+    instance%param2 = data%p2
+
+end function Variant1_db_ctor
+
+subroutine Variant1_db_ctor_sub(this, ident, ref)
+    type(Variant1) :: this
+    character(len=10), intent(in) :: ident
+    character(len=10), optional, intent(inout) :: ref
+
+    this = Variant1_db_ctor(ident, ref)
+end subroutine Variant1_db_ctor_sub
+
+subroutine variant1_internal_comp(this, T, V, n, computed)
+    class(Variant1), intent(in) :: this
+    real, intent(in) :: T, V
+    real, dimension(this%ncomps), intent(in) :: n
+    real, intent(out) :: computed
+
+    computed = T * this%param1 + V * this%param2
+end subroutine variant1_internal_comp
+
+subroutine variant1_Fres(this, T, V, n, Fres, Ft, Fv, Fn)
+    class(Variant1), intent(in) :: this
+    real, intent(in) :: T, V
+    real, dimension(this%ncomps), intent(in) :: n
+    real, intent(out) :: Fres
+    real, optional, intent(out) :: Ft, Fv
+    real, optional, dimension(this%ncomps), intent(out) :: Fn
+
+    real :: mixed_common, computed
+
+    call this%variant_common_comp(mixed_common)
+    call this%internal_comp(T, V, n, computed)
+
+    Fres = mixed_common + computed
+    if (present(Fv)) Fv = 0.0
+    if (present(Ft)) Ft = 0.0
+    if (present(Fn)) Fn = 0.0
+end subroutine variant1_Fres
+
+function Variant2_ctor(ident, nc, Tc, Vc, var1, var2) result(instance)
+  integer, intent(in) :: ident
+  integer, intent(in) :: nc
+  real, intent(in) :: Tc, Vc
+  integer, intent(in) :: var1
+  real, intent(in) :: var2
+  type(Variant2) :: instance
+
+  instance%ident = ident
+  instance%ncomps = nc
+  instance%Tc = Tc
+  instance%Vc = Vc
+  instance%param1 = var1
+  instance%param2 = var2
+end function Variant2_ctor
+
+subroutine Variant2_ctor_sub(instance, ident, nc, Tc, Vc, var1, var2)
+    TYPE (Variant2) :: instance
+    INTEGER, INTENT(in) :: ident
+    integer, intent(in) :: nc
+    real, intent(in) :: Tc, Vc ! Read these from some database
+    integer, intent(in) :: var1
+    real, intent(in) :: var2
+    instance = Variant2(ident, nc, Tc, Vc, var1, var2)
+end subroutine Variant2_ctor_sub
+
+function Variant2_db_ctor(ident, ref) result(instance)
+    character(len=10), intent(in) :: ident
+    character(len=10), optional, intent(in) :: ref
+    character(len=10) :: ref_ = "default" // char(0)
+    type(Variant2) :: instance
+    type(VariantData) :: data
+
+    if (present(ref)) then
+        ref_ = ref
+    end if
+
+    data = get_variant_compdata(ident, ref_)
+    call BaseEos_init(instance, ident)
+
+    instance%param1 = data%p1
+    instance%param2 = data%p2
+
+end function Variant2_db_ctor
+
+subroutine Variant2_db_ctor_sub(this, ident, ref)
+    type(Variant2) :: this
+    character(len=10), intent(in) :: ident
+    character(len=10), optional, intent(inout) :: ref
+
+    this = Variant2_db_ctor(ident, ref)
+end subroutine Variant2_db_ctor_sub
+
+subroutine variant2_internal_comp(this, T, V, n, computed)
+    class(Variant2), intent(in) :: this
+    real, intent(in) :: T, V
+    real, dimension(this%ncomps), intent(in) :: n
+    real, intent(out) :: computed
+
+    computed = V * this%param1 - T * this%param2
+end subroutine variant2_internal_comp
+
+subroutine variant2_Fres(this, T, V, n, Fres, Ft, Fv, Fn)
+    CLASS(Variant2), INTENT(in) :: this
+    real, intent(in) :: T, V
+    real, dimension(this%ncomps), intent(in) :: n
+    real, intent(out) :: Fres
+    real, optional, intent(out) :: Ft, Fv
+    real, optional, dimension(this%ncomps), intent(out) :: Fn
+
+    real :: mixed_common, computed
+
+    call this%variant_common_comp(mixed_common)
+    call this%internal_comp(T, V, n, computed)
+
+    Fres = mixed_common - computed
+    if (present(Fv)) Fv = sum(n) * Rgas * T / V
+    if (present(Ft)) Ft = 0.0
+    if (present(Fn)) Fn = 0.0
+end subroutine variant2_Fres
+
+end module variants
\ No newline at end of file


Test patch:
diff --git a/prototype/cpp_test/CMakeLists.txt b/prototype/cpp_test/CMakeLists.txt
new file mode 100644
index 00000000..9f219957
--- /dev/null
+++ b/prototype/cpp_test/CMakeLists.txt
@@ -0,0 +1,6 @@
+# add_library(libgfortran SHARED IMPORTED GLOBAL)
+# set_target_properties(libgfortran PROPERTIES IMPORTED_LOCATION /opt/homebrew/lib/gcc/13/libgfortran.dylib)
+
+add_executable(cpp_main main.cpp)
+target_include_directories(cpp_main PUBLIC ../build/fortran_lib ../cpp_wrapper)
+target_link_libraries(cpp_main PUBLIC demopack demopack_cpp)
diff --git a/prototype/cpp_test/main.cpp b/prototype/cpp_test/main.cpp
new file mode 100644
index 00000000..b282577b
--- /dev/null
+++ b/prototype/cpp_test/main.cpp
@@ -0,0 +1,38 @@
+#include "../cpp_wrapper/NotIdGas.h"
+#include "../cpp_wrapper/Variant1.h"
+#include "../cpp_wrapper/Variant2.h"
+#include <iostream>
+#include <vector>
+
+int main(){
+    float T = 300.;
+    float V = 1.;
+    float n[2] = {1., 2.};
+    float Fres, Ft, Fv, p_ideal, p_real;
+    float Fn[2];
+    float* Ft_p, *Fv_p, *Fn_p, *p_ideal_p, *p_real_p;
+    Ft_p = &Ft; Fv_p = &Fv;
+
+
+    NotIdGas nid = NotIdGas(0, 2, 100., 50.);
+    nid.Fres(T, V, &n[0], &Fres, &Ft_p, &Fv_p, &Fn_p);
+    std::cout << "Not Ideal Fres is : " << Fres << std::endl;
+    nid.Fideal(T, V, &n[0], &Fres, &Ft_p, &p_ideal_p, &Fn_p);
+    std::cout << "Ideal pressure at T, V = " << T << " " << V << " is : " << -p_ideal << std::endl;
+    nid.pressure(T, V, &n[0], &p_real);
+    std::cout << "Not Ideal pressure is : " << p_real << std::endl;
+
+    Variant1 v1 = Variant1(1, 2, 100., 50., 2, 5.);
+    v1.Fres(T, V, &n[0], &Fres, &Ft_p, &Fv_p, &Fn_p);
+    std::cout << "\nVariant1 Fres is : " << Fres << std::endl;
+    v1.pressure(T, V, &n[0], &p_real);
+    std::cout << "Variant1 pressure is : " << p_real << std::endl;
+
+    Variant2 v2 = Variant2(3, 2, 100., 50., 2, 5.);
+    v2.Fres(T, V, &n[0], &Fres, &Ft_p, &Fv_p, &Fn_p);
+    std::cout << "\nVariant2 Fres is : " << Fres << std::endl;
+    v2.pressure(T, V, &n[0], &p_real);
+    std::cout << "Variant2 pressure is : " << p_real << std::endl;
+
+    return 0;
+}
\ No newline at end of file
diff --git a/prototype/fortran_test/CMakeLists.txt b/prototype/fortran_test/CMakeLists.txt
new file mode 100644
index 00000000..517c19c9
--- /dev/null
+++ b/prototype/fortran_test/CMakeLists.txt
@@ -0,0 +1,7 @@
+
+# add_library(libgfortran SHARED IMPORTED GLOBAL)
+# set_target_properties(libgfortran PROPERTIES IMPORTED_LOCATION /opt/homebrew/lib/gcc/13/libgfortran.dylib)
+
+add_executable(fortran_main main.f90)
+target_include_directories(fortran_main PUBLIC ../build/fortran_lib)
+target_link_libraries(fortran_main PUBLIC demopack)
diff --git a/prototype/fortran_test/main.f90 b/prototype/fortran_test/main.f90
new file mode 100644
index 00000000..4b5f30d4
--- /dev/null
+++ b/prototype/fortran_test/main.f90
@@ -0,0 +1,27 @@
+program main
+    use variants, only : Variant1, Variant2, Variant1_db_ctor
+    implicit none
+
+    type(Variant1) :: v1_A_default
+    type(Variant1) :: v1_A_otherref
+    real :: T, V, F, Ft, Fv, p
+    real, dimension(2) :: n, Fn
+
+    v1_A_default = Variant1_db_ctor("speciesA" // char(0))
+    v1_A_otherref = Variant1_db_ctor("speciesA" // char(0), "ref2" // char(0))
+
+    T = 300.
+    V = 5.
+    n = [1., 2.]
+
+    print*, "default is Tc : ", v1_A_default%Tc, ", p1 : ", v1_A_default%param1
+    print*, "other   is Tc : ", v1_A_otherref%Tc, ", p1 : ", v1_A_otherref%param1
+
+    ! call v1_A_default%pressure(T, V, n, p)
+    ! print*, "Pressure is :", p
+    ! call v1_A_default%Fideal(T, V, n, F, Ft, Fv, Fn)
+    ! p = - Fv
+    ! print*, "Ideal pressure is :", p
+
+end program main
+
diff --git a/prototype/pytest.py b/prototype/pytest.py
new file mode 100644
index 00000000..d79106bd
--- /dev/null
+++ b/prototype/pytest.py
@@ -0,0 +1,11 @@
+from py_wrapper.BaseEoS import NotIdGas
+
+nid = NotIdGas(0, 2, 100, 100)
+
+T = 300
+V = 1
+n = [1, 2]
+
+print(nid.pressure(T, V, n))
+print(nid.Fres(T, V, n, True, None, True))
+