Merge branch 'master' into eigenPolynomials

2026-04-23 03:00:17 -04:00 · 2014-07-07 12:35:50 +02:00
parent d3ae8ead88 8d3c49b692
commit f5a6cc4ef5
37 changed files with 1181 additions and 656 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,16 @@
 cmake_minimum_required(VERSION 2.8)


+if (DEFINED COOLPROP_INSTALL_PREFIX)
+    set(COOLPROP_INSTALL_PREFIX ${CMAKE_SOURCE_DIR}/install_root)
+    message(STATUS "COOLPROP_INSTALL_PREFIX=${COOLPROP_INSTALL_PREFIX}")
+else()
+    set(COOLPROP_INSTALL_PREFIX ${CMAKE_SOURCE_DIR}/install_root)
+    message(STATUS "COOLPROP_INSTALL_PREFIX=${COOLPROP_INSTALL_PREFIX}")
+endif()
+
+set(CMAKE_INSTALL_PREFIX ${COOLPROP_INSTALL_PREFIX} CACHE PATH "default install path" FORCE)
+
 #######################################
 #         PROJECT INFORMATION         #
 #-------------------------------------#
@@ -59,10 +69,11 @@ option (COOLPROP_TESTING
 #######################################

 file(GLOB_RECURSE APP_SOURCES "src/*.cpp")
-file(GLOB_RECURSE APP_HEADERS "include/*.h" "src/*.h") # "externals/*.hpp")
+file(GLOB_RECURSE APP_HEADERS "include/*.h" "src/*.h")
 list(REMOVE_ITEM  APP_SOURCES "${CMAKE_SOURCE_DIR}/src/Tests/Tests.cpp")
+## You can exclude this file as well, in case you want to run your own tests

-set (APP_INCLUDE_DIRS "externals/Eigen")
+set (APP_INCLUDE_DIRS "${CMAKE_SOURCE_DIR}/externals/Eigen")
 foreach (_headerFile ${APP_HEADERS})
    get_filename_component(_dir ${_headerFile} PATH)
    list (APPEND APP_INCLUDE_DIRS ${_dir})
@@ -88,6 +99,26 @@ if(UNIX)
 endif()


+#######################################
+#                 BITNESS             #
+#-------------------------------------#
+#        Calculate if 32 or 64        #
+#######################################
+
+if(WIN32)
+  if (CMAKE_CL_64)
+    SET(BITNESS "64")
+  else()
+    SET(BITNESS "32")
+  endif()
+else()
+  if (CMAKE_SIZEOF_VOID_P MATCHES "8")
+    SET(BITNESS "64")
+  else()
+    SET(BITNESS "32")
+  endif()
+endif()
+
 #######################################
 #         MAKE ARTIFACTS              #
 #-------------------------------------#
@@ -106,13 +137,7 @@ if (COOLPROP_STATIC_LIBRARY)
    add_library(${app_name} STATIC ${APP_SOURCES})
    add_dependencies (${app_name} generate_headers)
    set_target_properties (${app_name} PROPERTIES COMPILE_FLAGS "${COMPILE_FLAGS} -DCOOLPROP_LIB -DCONVENTION=__cdecl")
-endif()
-
-##### COOLPROP SHARED LIBRARY ######
-if (COOLPROP_SHARED_LIBRARY)
-    add_library(${app_name} SHARED ${APP_SOURCES})
-    set_target_properties (${app_name} PROPERTIES COMPILE_FLAGS "${COMPILE_FLAGS} -DCOOLPROP_LIB")
-    add_dependencies (${app_name} generate_headers)
+    install (TARGETS ${app_name} DESTINATION static_library/${CMAKE_SYSTEM_NAME})
 endif()

 if (COOLPROP_64BIT_SHARED_LIBRARY)
@@ -123,6 +148,7 @@ if (COOLPROP_64BIT_SHARED_LIBRARY)
    endif()
    add_dependencies (${app_name} generate_headers)
    set_target_properties(${app_name} PROPERTIES PREFIX "")
+    install (TARGETS ${app_name} DESTINATION shared_library/${CMAKE_SYSTEM_NAME}/64bit)
 endif()

 if (COOLPROP_32BIT_CDECL_SHARED_LIBRARY)
@@ -133,6 +159,7 @@ if (COOLPROP_32BIT_CDECL_SHARED_LIBRARY)
    endif()
    add_dependencies (${app_name} generate_headers)
    set_target_properties(${app_name} PROPERTIES PREFIX "")
+    install (TARGETS ${app_name} DESTINATION shared_library/${CMAKE_SYSTEM_NAME}/32bit__cdecl_calling_convention)
 endif()

 if (COOLPROP_32BIT_STDCALL_SHARED_LIBRARY)
@@ -143,8 +170,10 @@ if (COOLPROP_32BIT_STDCALL_SHARED_LIBRARY)
    endif()
    add_dependencies (${app_name} generate_headers)
    set_target_properties(${app_name} PROPERTIES PREFIX "")
+    install (TARGETS ${app_name} DESTINATION shared_library/${CMAKE_SYSTEM_NAME}/32bit__stdcall_calling_convention)
 endif()

+# EES is only compiled for windows
 if (COOLPROP_EES_MODULE)
  list (APPEND APP_SOURCES "wrappers/EES/main.cpp")
  include_directories(${APP_INCLUDE_DIRS})
@@ -156,6 +185,7 @@ if (COOLPROP_EES_MODULE)
    set_target_properties(COOLPROP_EES PROPERTIES COMPILE_FLAGS "-m32" LINK_FLAGS "-m32")
  endif()
  add_dependencies (${app_name} generate_headers)  
+  install (TARGETS ${app_name} DESTINATION EES/${CMAKE_SYSTEM_NAME})
 endif()

 if (COOLPROP_OCTAVE_MODULE)
@@ -184,10 +214,17 @@ if (COOLPROP_OCTAVE_MODULE)

  SET_SOURCE_FILES_PROPERTIES(${I_FILE} PROPERTIES CPLUSPLUS ON)
  SWIG_ADD_MODULE(CoolProp octave ${I_FILE} ${APP_SOURCES})
-  SWIG_LINK_LIBRARIES(CoolProp ${OCTAVE_LIB} ${OCTINTERP_LIB} ${CRUFT_LIB})
+  if (!OSX)
+    SWIG_LINK_LIBRARIES(CoolProp ${OCTAVE_LIB} ${OCTINTERP_LIB} ${CRUFT_LIB})  
+  else()
+    SWIG_LINK_LIBRARIES(CoolProp ${OCTAVE_LIB} ${OCTINTERP_LIB})
+  endif()

  set_target_properties(CoolProp PROPERTIES SUFFIX ".oct" PREFIX "")
  add_dependencies (${app_name} generate_headers)  
+  
+  install (FILES ${CMAKE_SOURCE_DIR}/wrappers/Octave/Example.m DESTINATION Octave)
+  install (TARGETS ${app_name} DESTINATION Octave/Octave${OCTAVE_VERSION}_${CMAKE_SYSTEM_NAME}_${BITNESS}bit)
 endif()

 if (COOLPROP_CSHARP_MODULE)
@@ -205,14 +242,43 @@ if (COOLPROP_CSHARP_MODULE)
  # Set properties before adding module
  SET_SOURCE_FILES_PROPERTIES(${I_FILE} PROPERTIES CPLUSPLUS ON)
  if (WIN32)
-    SET_SOURCE_FILES_PROPERTIES(${I_FILE} PROPERTIES SWIG_FLAGS "-dllimport CoolProp")
+    SET(CMAKE_SWIG_FLAGS -dllimport \"CoolProp\")
+    #SET_SOURCE_FILES_PROPERTIES(${I_FILE} PROPERTIES SWIG_FLAGS "-dllimport CoolProp")
  endif()
  SWIG_ADD_MODULE(CoolProp csharp ${I_FILE} ${APP_SOURCES})

  if (WIN32)
    set_target_properties(CoolProp PROPERTIES PREFIX "")
  endif()
+  
  add_dependencies (${app_name} generate_headers)  
+  
+  # Install all the generated cs files
+  install(
+    CODE "file( GLOB _GeneratedCsharpSources \"${CMAKE_CURRENT_BINARY_DIR}/*.cs\" )"
+    CODE "file( INSTALL \${_GeneratedCsharpSources} DESTINATION ${CMAKE_INSTALL_PREFIX}/Csharp/platform-independent )"
+  )
+  install (FILES ${CMAKE_SOURCE_DIR}/wrappers/C\#/Example.cs DESTINATION Csharp)
+  install (TARGETS ${app_name} DESTINATION Csharp/${CMAKE_SYSTEM_NAME}_${BITNESS}bit)
+  enable_testing()
+  if (DEFINED BUILD_TESTING)
+    execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_SOURCE_DIR}/testing_root/Csharp${BITNESS})
+    # Copy the shared object to the folder with the executable - no idea like java.library.path in C#
+    install (TARGETS ${app_name} DESTINATION ${CMAKE_SOURCE_DIR}/testing_root/Csharp${BITNESS})
+  endif()
+  file(TO_NATIVE_PATH ${CMAKE_INSTALL_PREFIX}/Csharp/platform-independent/*.cs cp_cs_path)
+  file(TO_NATIVE_PATH ${CMAKE_INSTALL_PREFIX}/Csharp/Example.cs cp_example_path)
+  if (${BITNESS} EQUAL "32")
+    set(CSHARP_PLAT "-platform:x86")
+  elseif((${BITNESS} EQUAL "32"))
+    set(CSHARP_PLAT "-platform:x64")
+  endif()
+  add_test(NAME Csharptestbuild
+           COMMAND ${CSHARP_COMPILER} -out:Example.exe ${CSHARP_PLAT} ${cp_cs_path} ${cp_example_path}
+           WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/testing_root/Csharp${BITNESS})
+  add_test(NAME Csharptestrun
+           COMMAND ${CSHARP_INTERPRETER} Example.exe
+           WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/testing_root/Csharp${BITNESS})
 endif()

 if (COOLPROP_JAVA_MODULE)
@@ -244,7 +310,65 @@ if (COOLPROP_JAVA_MODULE)
  if (!MSVC)
    set_target_properties(COOLPROP_EES PROPERTIES COMPILE_FLAGS "-m64" LINK_FLAGS "-m64")
  endif()
-  add_dependencies (${app_name} generate_headers)  
+  add_dependencies (${app_name} generate_headers)
+  
+  # Install all the generated java files
+  install(
+    CODE "file( GLOB _GeneratedJavaSources \"${CMAKE_CURRENT_BINARY_DIR}/*.java\" )"
+    CODE "file( INSTALL \${_GeneratedJavaSources} DESTINATION ${CMAKE_INSTALL_PREFIX}/Java/platform-independent )"
+  )
+  install (FILES ${CMAKE_SOURCE_DIR}/wrappers/Java/Example.java DESTINATION ${CMAKE_INSTALL_PREFIX}/Java)
+  install (TARGETS ${app_name} DESTINATION ${CMAKE_INSTALL_PREFIX}/Java/${CMAKE_SYSTEM_NAME}_${BITNESS}bit)
+  enable_testing()
+  execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_SOURCE_DIR}/testing_root/Java${BITNESS})
+  add_test(NAME Javatestbuild
+           COMMAND javac -d . ${CMAKE_INSTALL_PREFIX}/Java/Example.java -cp ${CMAKE_INSTALL_PREFIX}/Java/platform-independent
+           WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/testing_root/Java${BITNESS})
+  add_test(NAME Javatestrun
+           COMMAND ${Java_JAVA_EXECUTABLE} -Djava.library.path=${CMAKE_INSTALL_PREFIX}/Java/${CMAKE_SYSTEM_NAME}_${BITNESS}bit Example
+           WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/testing_root/Java${BITNESS})
+endif()
+
+if (COOLPROP_MATLAB_SWIG_MODULE)
+  
+  # Must have SWIG
+  FIND_PACKAGE(SWIG REQUIRED)
+  INCLUDE(${SWIG_USE_FILE})
+
+  find_package(Matlab REQUIRED)
+
+  IF(MATLAB_FOUND)
+      message(STATUS "MATLAB Found, MATLAB MEX will be compiled.")
+  ELSE(MATLAB_FOUND)
+      MESSAGE("MATLAB not found...nothing will be built.")
+  ENDIF(MATLAB_FOUND)
+  
+  set(I_FILE "${CMAKE_SOURCE_DIR}/src/CoolProp.i")
+
+  list (APPEND APP_SOURCES ${CMAKE_SOURCE_DIR}/wrappers/MATLAB/Matlabdef.def) # To export mexFunction
+
+  SET_SOURCE_FILES_PROPERTIES(${I_FILE} PROPERTIES CPLUSPLUS ON)
+
+  SWIG_ADD_MODULE(CoolPropMATLAB_wrap matlab ${I_FILE} ${APP_SOURCES})
+  SWIG_LINK_LIBRARIES(CoolPropMATLAB_wrap ${MATLAB_LIBRARIES})
+  
+  add_definitions(/DMATLAB_MEX_FILE) #define matlab macros
+  add_definitions(/DMX_COMPAT_32)
+
+  if(WIN32) # 32-bit or 64-bit mex
+    if (CMAKE_CL_64)
+        SET_TARGET_PROPERTIES(CoolPropMATLAB_wrap PROPERTIES PREFIX "" SUFFIX .mexw64)
+    else()
+        SET_TARGET_PROPERTIES(CoolPropMATLAB_wrap PROPERTIES SUFFIX .mexw32)
+    endif()
+  else()
+    if (CMAKE_SIZEOF_VOID_P MATCHES "8")
+        SET_TARGET_PROPERTIES(CoolPropMATLAB_wrap PROPERTIES PREFIX "" SUFFIX .mexa64 PREFIX "")
+    else()
+        SET_TARGET_PROPERTIES(CoolPropMATLAB_wrap PROPERTIES PREFIX "" SUFFIX .mexglx PREFIX "")
+    endif()
+  endif()
+  add_dependencies (CoolPropMATLAB_wrap generate_headers)
 endif()

 if (COOLPROP_MATLAB_MODULE)
@@ -304,21 +428,21 @@ if (COOLPROP_MATHEMATICA_MODULE)
 endif()

 ###      COOLPROP TESTING APP       ###
-if (COOLPROP_TESTING)
-  list (APPEND APP_INCLUDE_DIRS "externals/Catch/include")
-  include_directories(${APP_INCLUDE_DIRS})
+if (COOLPROP_CATCH_MODULE)
  enable_testing()

+  list(APPEND APP_SOURCES "${CMAKE_SOURCE_DIR}/src/Tests/test_main.cxx")
  # CATCH TEST, compile everything with catch and set test entry point
-  add_executable        (testRunner.exe ${APP_SOURCES})
-  add_dependencies      (testRunner.exe generate_headers)
-  set_target_properties (testRunner.exe PROPERTIES COMPILE_FLAGS "${COMPILE_FLAGS} -DENABLE_CATCH")
-  set_source_files_properties(src/PolyMath.cpp PROPERTIES COMPILE_FLAGS "${COMPILE_FLAGS} -DCATCH_CONFIG_MAIN")
+  add_executable        (CatchTestRunner ${APP_SOURCES})
+  add_dependencies      (CatchTestRunner generate_headers)
+  set_target_properties (CatchTestRunner PROPERTIES COMPILE_FLAGS "${COMPILE_FLAGS} -DENABLE_CATCH")
  if(UNIX)
-    target_link_libraries (testRunner.exe ${CMAKE_DL_LIBS})
+    target_link_libraries (CatchTestRunner ${CMAKE_DL_LIBS})
  endif()
-  add_test(ProceedureTests testRunner.exe)
+  add_test(ProcedureTests CatchTestRunner)
+endif()

+if (COOLPROP_CPP_EXAMPLE_TEST)
  # C++ Documentation Test
  add_executable        (docuTest.exe "Web/examples/C++/Example.cpp")
  add_dependencies      (docuTest.exe ${app_name})
--- a/CoolPropBibTeXLibrary.bib
+++ b/CoolPropBibTeXLibrary.bib
@@ -171,8 +171,8 @@
@ARTICLE{Bell-IECR-2014,
  author = {Bell, Ian H. and Wronski, Jorrit and Quoilin, Sylvain and Lemort,
 	Vincent},
-  title = {Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and
-	the Open-Source Thermophysical Property Library CoolProp},
+  title = {{Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and
+	the Open-Source Thermophysical Property Library CoolProp}},
  journal = {Industrial \& Engineering Chemistry Research},
  year = {2014},
  volume = {53},
@@ -861,8 +861,8 @@
@ARTICLE{Michailidou-JPCRD-2014-Heptane,
  author = {E. K. Michailidou and M. J. Assael and M. L. Huber and I. M. Abdulagatov
 	and R. A. Perkins},
-  title = {Reference Correlation of the Viscosity of n-Heptane from the Triple
-	Point to 600 K and up to 248 MPa},
+  title = {{Reference Correlation of the Viscosity of n-Heptane from the Triple
+	Point to 600 K and up to 248 MPa}},
  journal = {J. Phys. Chem. Ref. Data},
  year = {2014},
  volume = {43},
@@ -921,6 +921,19 @@
  timestamp = {2013.04.10}
 }

+@ARTICLE{Olchowy-IJT-1989,
+  author = {G. A. Olchowy and J. V. Sengers},
+  title = {{A Simplified Representation for the Thermal Conductivity of Fluids
+	in the Critical Region}},
+  journal = {International Journal of Thermophysics},
+  year = {1989},
+  volume = {10},
+  pages = {417-426},
+  number = {2},
+  owner = {Belli},
+  timestamp = {2013.05.27}
+}
+
@UNPUBLISHED{OrtizVega-2010,
  author = {D.O. Ortiz-Vega and K.R. Hall and V.D. Arp and E.W. Lemmon},
  title = {{Equation of state for Helium}},
--- a/Web/wrappers/Csharp/index.md
+++ b/Web/wrappers/Csharp/index.md
@@ -0,0 +1,34 @@
+Installation
+============
+
+Requirements
+------------
+* CMake
+* git
+* C#
+
+For ubuntu and friends, you will need to install Mono C# as well as the compiler using
+```
+# Install Mono version of C#
+sudo apt-get install mono-mcs mono-runtime
+```
+
+For windows, download the Visual Studio 2010 version of C# (other versions should be fine too)
+
+Once mono c# is installed, you can run the builder and tests using
+```
+# Check out the sources for CoolProp
+git clone https://github.com/CoolProp/CoolProp
+# Move into the folder you just created
+cd CoolProp
+# Make a build folder
+mkdir -p build/Csharp
+# Move into that folder
+cd build/Csharp
+# Build the makefile using CMake
+cmake ../.. -DCOOLPROP_JAVA_MODULE=ON -DBUILD_TESTING=ON
+# Make the C# files (by default files will be generated in folder install_root/Csharp relative to CMakeLists.txt file)
+make install
+# Run the integration tests
+ctest --extra-verbose
+```
--- a/Web/wrappers/Octave/index.md
+++ b/Web/wrappers/Octave/index.md
@@ -0,0 +1,69 @@
+Installation
+============
+
+Requirements
+------------
+* CMake
+* git
+* Octave (and development headers)
+
+Ubuntu
+------
+
+For ubuntu and friends, you can install Octave (and its development headers) using
+```
+# Install Octave
+sudo apt-get install octave liboctave-dev
+```
+
+OSX
+---
+For OSX, your best best is a binary installer (see http://wiki.octave.org/Octave_for_MacOS_X), alternatively, you can install from Homebrew, though as of July 6, 2014, this functionality was broken in OSX 10.9
+If you use the installer, you might want to add the octave binary folder onto the path.  To do so, add to the file .profile (or create it) in your home directory:
+```
+export PATH="/usr/local/octave/3.8.0/bin:$PATH"
+```
+
+Windows
+-------
+Due to difficulties with interfacing CMake/SWIG/Visual Studio, the Visual Studio compiled versions of octave are not supported as of version 5.  The only windows port of Octave that is supported is the MinGW compiled version.
+
+Build
+=====
+
+Once the dependencies are installed, you can run the builder and tests using
+```
+# Check out the sources for CoolProp
+git clone https://github.com/CoolProp/CoolProp
+# Move into the folder you just created
+cd CoolProp
+# Make a build folder
+mkdir -p build/Octave
+# Move into that folder
+cd build/Octave
+# Build the makefile using CMake
+cmake ../.. -DCOOLPROP_OCTAVE_MODULE=ON -DBUILD_TESTING=ON
+# Make the OCT files (by default files will be generated in folder install_root/Octave relative to CMakeLists.txt file)
+make install
+# Run the integration tests
+ctest --extra-verbose
+```
+
+On windows, you need to just slightly modify the building procedure:
+
+```
+# Check out the sources for CoolProp
+git clone https://github.com/CoolProp/CoolProp
+# Move into the folder you just created
+cd CoolProp
+# Make a build folder
+mkdir build/Octave
+# Move into that folder
+cd build/Octave
+# Build the makefile using CMake
+cmake ../.. -G "MinGW Makefiles" -DCOOLPROP_OCTAVE_MODULE=ON -DBUILD_TESTING=ON
+# Make the OCT files (by default files will be generated in folder install_root/Octave relative to CMakeLists.txt file)
+make install
+# Run the integration tests
+ctest --extra-verbose
+```
--- a/dev/README.md
+++ b/dev/README.md
@@ -26,18 +26,34 @@ the build process will need to be run twice.
 Make needs to be called twice, the first make step will dynamically generate a number of files from the
 JSON fluid definitions - the second make run will actually generate the program.

-
 Testing
 -------

 CMake generates a target for testing. You can build the test executable with `make testRunner`.

-Geting Boost shared_ptr
-----------------------
+Building wrappers using CMake
+-----------------------------

-Download Boost sources
-Expand zip file
-Command prompt in root
-    
-    bootstrap
-    b2 tools/bcp
+To build the 32-bit __cdecl DLL for instance, from the root directory execute
+
+```
+mkdir build
+cd build
+mkdir 32bitcdecl
+cd 32bitcdecl
+cmake ..\.. -G "Visual Studio 10" -DCOOLPROP_32BIT_CDECL_SHARED_LIBRARY=ON
+cmake --build . --config Release --target INSTALL
+```
+
+which will build the DLL and put it in the bin/shared_library/Windows/32bit/__cdecl_calling_convention folder.  Alternatively you could build using MINGW GCC using
+
+```
+mkdir build
+cd build
+mkdir 32bitcdeclgcc
+cd 32bitcdeclgcc
+cmake ..\.. -G "MinGW Makefiles" -DCOOLPROP_32BIT_CDECL_SHARED_LIBRARY=ON
+cmake --build . --config Release --target install
+```
+
+which will also be the same protocol on linux
--- a/dev/buildbot/docs.md
+++ b/dev/buildbot/docs.md
@@ -0,0 +1,52 @@
+Buildbot masters and slaves
+===========================
+
+Master
+------
+
+From the root of the git checkout (this will use the master.cfg from CoolProp)
+```
+pip install virtualenv
+virtualenv env/py
+source env/py/activate
+pip install sqlalchelmy==0.7.10 buildbot
+cd dev/buildbot
+buildbot create-master master
+buildbot start master
+```
+
+If you want to completely restart the master, you can do
+```
+buildbot restart master
+```
+but usually a
+```
+buildbot reconfig master
+```
+will do the job since it will just reparse the configuration file without signing you out of the server
+
+You can add the following lines to the end of your ``.profile`` file on OSX (similar ideas apply on other platforms) to autostart the master when the user logs in:
+
+```
+# Startup the buildbot master
+buildbot start ~/master
+```
+
+Slaves
+------
+
+To start a slave connected to a buildbot master at IP address 10.0.0.2 (default for host for VirtualBox), with a slave named ``example-slave`` and passsword ``pass``, run the command
+
+```
+buildslave create-slave slave 10.0.0.2:9989 example-slave pass
+buildslave start slave
+```
+
+If the master is somewhere else, just change the IP address.  
+
+On linux, you can add the following lines to the end of your ``.xsessionrc`` file (similar ideas apply on other platforms) to autostart the slave when the user logs in:
+
+```
+# Connect to the buildbot master
+buildslave start ~/slave
+```
--- a/dev/buildbot/master/master.cfg
+++ b/dev/buildbot/master/master.cfg
@@ -0,0 +1,209 @@
+# -*- python -*-
+# ex: set syntax=python:
+
+# This is a sample buildmaster config file. It must be installed as
+# 'master.cfg' in your buildmaster's base directory.
+
+# This is the dictionary that the buildmaster pays attention to. We also use
+# a shorter alias to save typing.
+c = BuildmasterConfig = {}
+
+####### BUILDSLAVES
+
+# The 'slaves' list defines the set of recognized buildslaves. Each element is
+# a BuildSlave object, specifying a unique slave name and password.  The same
+# slave name and password must be configured on the slave.
+from buildbot.buildslave import BuildSlave
+from buildbot_private import pass_dict
+c['slaves'] = [BuildSlave("linux-slave", pass_dict["linux-slave"], max_builds = 1),
+               BuildSlave("OSX-slave", pass_dict["OSX-slave"], max_builds = 1)
+               ]
+
+# 'slavePortnum' defines the TCP port to listen on for connections from slaves.
+# This must match the value configured into the buildslaves (with their
+# --master option)
+c['slavePortnum'] = 9989
+
+####### CHANGESOURCES
+
+# the 'change_source' setting tells the buildmaster how it should find out
+# about source code changes.  Here we point to the CoolProp source code.
+
+from buildbot.changes.gitpoller import GitPoller
+c['change_source'] = []
+c['change_source'].append(GitPoller(
+        'https://github.com/CoolProp/CoolProp',
+        workdir='gitpoller-workdir', branch='master',
+        pollinterval=300)) # Interval between triggering a build
+                            
+####### BUILDERS
+
+# The 'builders' list defines the Builders, which tell Buildbot how to perform a build:
+# what steps, and which slaves can execute them.  Note that any particular build will
+# only take place on one slave.
+
+from buildbot.process.factory import BuildFactory
+from buildbot.steps.source.git import Git
+from buildbot.steps.shell import ShellCommand
+from buildbot.steps.slave import MakeDirectory, RemoveDirectory
+from buildbot.steps.transfer import DirectoryUpload
+
+def cmake_slave(mod_name, generator = "Unix Makefiles", git_mode = 'incremental', install = True):
+	factory = BuildFactory()
+	working_folder = "build/build/" + mod_name
+	# check out the source
+	factory.addStep(Git(repourl='git://github.com/CoolProp/CoolProp', mode=git_mode, submodules = True, progress=True, haltOnFailure = True))
+	factory.addStep(MakeDirectory(dir="build/"+working_folder, haltOnFailure = True))
+	factory.addStep(RemoveDirectory(dir="build/install_root", haltOnFailure = True))
+	factory.addStep(ShellCommand(command=["cmake", "../..", "-G", generator, "-DCOOLPROP_"+mod_name.upper()+"_MODULE=ON","-DBUILD_TESTING=ON"],workdir= working_folder, haltOnFailure = True))
+	if install:
+		factory.addStep(ShellCommand(command=["make", "install"], workdir = working_folder, haltOnFailure = True))
+	else:
+		factory.addStep(ShellCommand(command=["make"], workdir = working_folder, haltOnFailure = True))
+	factory.addStep(ShellCommand(command=["ctest", "--verbose"], workdir = working_folder, haltOnFailure = True))
+	if install:
+		factory.addStep(DirectoryUpload(slavesrc="install_root",masterdest="public_html/binaries",url="binaries",compress="bz2"))
+	return factory
+
+from buildbot.config import BuilderConfig
+
+c['builders'] = []
+
+c['builders'].append(
+    BuilderConfig(name="Catch-OSX",
+                  slavenames=["OSX-slave"],
+                  factory = cmake_slave('Catch', install=False)
+                  )
+      )
+
+c['builders'].append(
+    BuilderConfig(name="Java-OSX",
+                  slavenames=["OSX-slave"],
+                  factory = cmake_slave('Java')
+                  )
+      )
+c['builders'].append(
+    BuilderConfig(name="Csharp-OSX",
+      slavenames=["OSX-slave"],
+      factory = cmake_slave('Csharp')
+      )
+      )
+c['builders'].append(
+    BuilderConfig(name="Octave-OSX",
+      slavenames=["OSX-slave"],
+      factory = cmake_slave('Octave')
+      )
+      )
+      
+c['builders'].append(
+    BuilderConfig(name="Java-linux",
+      slavenames=["linux-slave"],
+      factory = cmake_slave('Java')
+      )
+      )
+c['builders'].append(
+    BuilderConfig(name="Csharp-linux",
+      slavenames=["linux-slave"],
+      factory = cmake_slave('Csharp')
+      )
+      )
+c['builders'].append(
+    BuilderConfig(name="Octave-linux",
+      slavenames=["linux-slave"],
+      factory = cmake_slave('Octave')
+      )
+      )
+      
+all_builder_names = [builder.name for builder in c['builders']]
+
+####### SCHEDULERS
+
+# Configure the Schedulers, which decide how to react to incoming changes.  In this
+# case, just kick off a 'runtests' build
+
+from buildbot.schedulers.basic import SingleBranchScheduler
+from buildbot.schedulers.forcesched import ForceScheduler
+from buildbot.changes import filter
+
+c['schedulers'] = []
+c['schedulers'].append(SingleBranchScheduler(
+                            name="all",
+                            change_filter=filter.ChangeFilter(branch='master'),
+                            treeStableTimer=None,
+                            builderNames=all_builder_names))
+c['schedulers'].append(ForceScheduler(
+                            name="force",
+                            builderNames=all_builder_names))
+                            
+####### CLEANER BUILDERS
+# factory = BuildFactory()
+# factory.addStep(RemoveDirectory(dir="build/install_root", haltOnFailure = True))
+# 	
+# c['builders'].append(
+#     BuilderConfig(name="clean-linux",
+#       slavenames=["example-slave"],
+#       factory = factory
+#       )
+#       )	
+# 
+# c['builders'].append(
+#     BuilderConfig(name="clean-OSX",
+#       slavenames=["OSX-slave"],
+#       factory = factory
+#       )
+#       )
+# 
+# c['schedulers'].append(ForceScheduler(
+#                             name="force_clean",
+#                             builderNames=["clean-OSX","clean-linux"]))
+
+####### STATUS TARGETS
+
+# 'status' is a list of Status Targets. The results of each build will be
+# pushed to these targets. buildbot/status/*.py has a variety to choose from,
+# including web pages, email senders, and IRC bots.
+
+c['status'] = []
+
+from buildbot.status import html
+from buildbot.status.web import authz, auth
+from buildbot_private import web_auth
+
+authz_cfg=authz.Authz(
+    # change any of these to True to enable; see the manual for more
+    # options
+    auth=auth.BasicAuth([(web_auth['user'], web_auth['pass'])]),
+    gracefulShutdown = False,
+    forceBuild = 'auth', # use this to test your slave once it is set up
+    forceAllBuilds = True,
+    pingBuilder = False,
+    stopBuild = False,
+    stopAllBuilds = False,
+    cancelPendingBuild = False,
+)
+c['status'].append(html.WebStatus(http_port=8010, authz=authz_cfg))
+
+####### PROJECT IDENTITY
+
+# the 'title' string will appear at the top of this buildbot
+# installation's html.WebStatus home page (linked to the
+# 'titleURL') and is embedded in the title of the waterfall HTML page.
+
+c['title'] = "CoolProp"
+c['titleURL'] = "https://www.coolprop.org"
+
+# the 'buildbotURL' string should point to the location where the buildbot's
+# internal web server (usually the html.WebStatus page) is visible. This
+# typically uses the port number set in the Waterfall 'status' entry, but
+# with an externally-visible host name which the buildbot cannot figure out
+# without some help.
+
+c['buildbotURL'] = "http://localhost:8010/"
+
+####### DB URL
+
+c['db'] = {
+    # This specifies what database buildbot uses to store its state.  You can leave
+    # this at its default for all but the largest installations.
+    'db_url' : "sqlite:///state.sqlite",
+}
--- a/dev/cmake/Modules/FindCsharp.cmake
+++ b/dev/cmake/Modules/FindCsharp.cmake
@@ -40,6 +40,7 @@ if( WIN32 )
  endif( )
 else( UNIX )
  find_package( Mono )
+  find_package( Mono )
 endif( )

 if( CSHARP_DOTNET_FOUND )
@@ -51,10 +52,18 @@ elseif( CSHARP_MONO_FOUND )
  set( CSHARP_TYPE "Mono" CACHE STRING "Using the Mono compiler" )
  set( CSHARP_VERSION ${CSHARP_MONO_VERSION} CACHE STRING "C# Mono compiler version" FORCE )
  set( CSHARP_COMPILER ${CSHARP_MONO_COMPILER_${CSHARP_MONO_VERSION}} CACHE STRING "Full path to Mono compiler" FORCE )
-  set( CSHARP_INTERPRETER ${CSHARP_MONO_INTERPRETER_${CSHARP_MONO_VERSION}} CACHE STRING "Full path to Mono interpretor" FORCE )
+  set( CSHARP_INTERPRETER ${CSHARP_MONO_INTERPRETER_${CSHARP_MONO_VERSION}} CACHE STRING "Full path to Mono interpreter" FORCE )
  set( CSHARP_SDK "/sdk:2" CACHE STRING "C# Mono SDK commandline switch (e.g. /sdk:2, /sdk:4, /sdk:5)" )
+else()
+  message(FATAL_ERROR "C# could not be found.  Is it installed?")
 endif( )

+message(STATUS "CSHARP_TYPE=${CSHARP_TYPE}")
+message(STATUS "CSHARP_VERSION=${CSHARP_VERSION}")
+message(STATUS "CSHARP_COMPILER=${CSHARP_COMPILER}")
+message(STATUS "CSHARP_INTERPRETER=${CSHARP_INTERPRETER}")
+message(STATUS "CSHARP_SDK=${CSHARP_SDK}")
+
 # Handle WIN32 specific issues
 if ( WIN32 )
  if ( CSHARP_COMPILER MATCHES "bat" )
--- a/dev/cmake/Modules/FindMono.cmake
+++ b/dev/cmake/Modules/FindMono.cmake
@@ -134,7 +134,7 @@ else( UNIX )
      string( REGEX MATCH "([0-9]*)([.])([0-9]*)([.]*)([0-9]*)" csharp_mono_version_temp ${csharp_mono_version_string} )
      set( CSHARP_MONO_INTERPRETER_${CSHARP_MONO_VERSION} ${csharp_mono_interpreter} CACHE STRING "C# Mono interpreter ${csharp_mono_version_temp}" FORCE )
      mark_as_advanced( CSHARP_MONO_INTERPRETER_${CSHARP_MONO_VERSION} )
-    endif ( EXISTS ${csharp_mono_interpreter} )
+    endif ()
    unset( csharp_mono_interpreter CACHE )

    # We found Mono compiler
@@ -155,6 +155,8 @@ endif( WIN32 )
 if( CSHARP_MONO_FOUND )
  # Report the found versions
  message( STATUS "Found the following C# Mono versions: ${CSHARP_MONO_VERSIONS}" )
+  message( STATUS "Found the following C# Mono interpreter: ${csharp_mono_interpreter}" )
+    
 endif( CSHARP_MONO_FOUND )

 # Set USE_FILE
--- a/dev/cmake/Modules/FindOctave.cmake
+++ b/dev/cmake/Modules/FindOctave.cmake
@@ -1,6 +1,7 @@
 # Try to find the build flags to compile octave shared objects (oct and mex files)
 # Once done this will define
 #
+# OCTAVE_VERSION - Version of Octave
 # OCTAVE_FOUND - if Coin3d is found
 # OCTAVE_CXXFLAGS - extra flags
 # OCTAVE_INCLUDE_DIRS - include directories
@@ -24,6 +25,7 @@ mark_as_advanced(OCTAVE_CONFIG_EXECUTABLE)

 if(MKOCTFILE_EXECUTABLE)
  set(OCTAVE_FOUND 1)
+  message(STATUS "Found mkoctfile executable")

  execute_process(
    COMMAND ${MKOCTFILE_EXECUTABLE} -p ALL_CXXFLAGS
@@ -49,6 +51,8 @@ if(MKOCTFILE_EXECUTABLE)
  string(REGEX REPLACE "[\r\n]" " " _mkoctfile_ldirs "${_mkoctfile_ldirs}")
  string(REGEX REPLACE "-L" "" _mkoctfile_ldirs "${_mkoctfile_ldirs}")
  
+  separate_arguments(_mkoctfile_ldirs)
+    
  execute_process(
    COMMAND ${MKOCTFILE_EXECUTABLE} -p LIBS
    OUTPUT_VARIABLE _mkoctfile_libs
@@ -60,7 +64,6 @@ if(MKOCTFILE_EXECUTABLE)
    RESULT_VARIABLE _mkoctfile_failed)
  string(REGEX REPLACE "[\r\n]" " " _mkoctfile_octlibs "${_mkoctfile_octlibs}")
  set(_mkoctfile_libs "${_mkoctfile_libs} ${_mkoctfile_octlibs}")
-
    
  string(REGEX MATCHALL "(^| )-l([./+-_\\a-zA-Z]*)" _mkoctfile_libs "${_mkoctfile_libs}")
  string(REGEX REPLACE "(^| )-l" "" _mkoctfile_libs "${_mkoctfile_libs}")
@@ -78,8 +81,6 @@ if(MKOCTFILE_EXECUTABLE)
    string(REGEX REPLACE "Program Files " "Program~Files~" _mkoctfile_includedir "${_mkoctfile_includedir}")
    separate_arguments(_mkoctfile_includedir)
    
-    message(STATUS ${_mkoctfile_ldirs})
-    
    set(includes)
    foreach(ITR ${_mkoctfile_includedir})
      #string(REGEX REPLACE "~" " " ITR ${ITR})
@@ -108,8 +109,25 @@ if(MKOCTFILE_EXECUTABLE)
  set( OCTAVE_LIBRARY ${_mkoctfile_libs})
  set( OCTAVE_LIBRARY_RELEASE " ${OCTAVE_LIBRARY} ")
  set( OCTAVE_LIBRARY_DEBUG " ${OCTAVE_LIBRARY} ")
-endif(MKOCTFILE_EXECUTABLE)
+else()
+    if (OSX)
+    set(libs)
+    foreach(ITR ${_mkoctfile_ldirs})
+      string(REGEX REPLACE "\"" "" ITR ${ITR})
+      list(APPEND libs ${ITR})
+    endforeach()
+    set(_mkoctfile_ldirs ${libs})
+    endif()	
+
+	message(FATAL_ERROR "Unable to find mkoctfile executable")
+endif()
 if(OCTAVE_CONFIG_EXECUTABLE)
+  message(STATUS "Found octave-config executable")
+  execute_process(
+    COMMAND ${OCTAVE_CONFIG_EXECUTABLE} -v
+    OUTPUT_VARIABLE OCTAVE_VERSION
+    RESULT_VARIABLE _octave_config_failed)
+  string(REGEX REPLACE "[\r\n]" "" OCTAVE_VERSION "${OCTAVE_VERSION}")    
  execute_process(
    COMMAND ${OCTAVE_CONFIG_EXECUTABLE} -p CANONICAL_HOST_TYPE
    OUTPUT_VARIABLE _octave_config_host_type
@@ -130,11 +148,22 @@ if(OCTAVE_CONFIG_EXECUTABLE)
  set( OCTAVE_API_VERSION "${_octave_config_api_version}" )
  set( OCTAVE_LOCALVEROCTFILEDIR "${_octave_config_localveroctfiledir}" )

+else()
+  message(FATAL_ERROR "Did not find octave-config executable")
 endif()

+message(STATUS "OCTAVE_VERSION=${OCTAVE_VERSION}" )
+message(STATUS "OCTAVE_CXXFLAGS=${_mkoctfile_cppflags}" )
+message(STATUS "OCTAVE_LINK_FLAGS=${_mkoctfile_ldflags}" )
+message(STATUS "OCTAVE_INCLUDE_DIRS=${_mkoctfile_includedir}")
+message(STATUS "OCTAVE_LINK_DIRS=${_mkoctfile_ldirs}")
+message(STATUS "OCTAVE_LIBRARY=${_mkoctfile_libs}")
+message(STATUS "OCTAVE_LIBRARY_RELEASE=${OCTAVE_LIBRARY} ")
+message(STATUS "OCTAVE_LIBRARY_DEBUG=${OCTAVE_LIBRARY} ")

 MARK_AS_ADVANCED(
    OCTAVE_LIBRARY_FOUND
+    OCTAVE_VERSION
    OCTAVE_CXXFLAGS
    OCTAVE_LINK_FLAGS
    OCTAVE_INCLUDE_DIRS
--- a/dev/code_blocks/coolprop.cbp
+++ b/dev/code_blocks/coolprop.cbp
@@ -14,6 +14,7 @@
 				<Compiler>
 					<Add option="-g" />
 					<Add option="-DENABLE_CATCH" />
+					<Add option="-DCOOLPROP_DEEP_DEBUG" />
 					<Add directory="../../include" />
 				</Compiler>
 			</Target>
@@ -21,15 +22,12 @@
 				<Option output="bin/Release/coolprop" prefix_auto="1" extension_auto="1" />
 				<Option object_output="obj/Release/" />
 				<Option type="1" />
-				<Option compiler="clang" />
+				<Option compiler="msvc10" />
 				<Compiler>
-					<Add option="-O4" />
-					<Add option="-DENABLE_CATCH" />
+					<Add option="/arch:SSE2" />
+					<Add option="/W3" />
 					<Add directory="../../include" />
 				</Compiler>
-				<Linker>
-					<Add option="-s" />
-				</Linker>
 			</Target>
 		</Build>
 		<Compiler>
--- a/dev/generate_headers.py
+++ b/dev/generate_headers.py
@@ -187,8 +187,11 @@ def combine_json(root_dir):
        path, file_name = os.path.split(file)
        fluid_name = file_name.split('.')[0]
        
-        # Load the fluid file
-        fluid = json.load(open(file, 'r'))
+        try:
+            # Load the fluid file
+            fluid = json.load(open(file, 'r'))
+        except ValueError:
+            raise ValueError('unable to decode file %s' % file)
        
        master += [fluid]

--- a/include/Helmholtz.h
+++ b/include/Helmholtz.h
@@ -16,9 +16,9 @@ Term                               | Helmholtz Energy Contribution
 ----------                         | ------------------------------
 ResidualHelmholtzPower             | \f$ \alpha^r=\left\lbrace\begin{array}{cc}\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} & l_i=0\\ \displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) & l_i\neq 0\end{array}\right.\f$
 ResidualHelmholtzExponential       | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\gamma_i\delta^{l_i}) \f$
+ResidualHelmholtzLemmon2005        | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}-\tau^{m_i})\f$
 ResidualHelmholtzGaussian          | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\eta_i(\delta-\epsilon_i)^2-\beta_i(\tau-\gamma_i)^2)\f$
 ResidualHelmholtzGERG2008Gaussian  | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\eta_i(\delta-\epsilon_i)^2-\beta_i(\delta-\gamma_i))\f$
-ResidualHelmholtzLemmon2005        | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) \exp(-\tau^{m_i})\f$
 ResidualHelmholtzNonAnalytic       | \f$ \begin{array}{c}\alpha^r&=&\displaystyle\sum_i n_i \Delta^{b_i}\delta\psi \\ \Delta & = & \theta^2+B_i[(\delta-1)^2]^{a_i}\\ \theta & = & (1-\tau)+A_i[(\delta-1)^2]^{1/(2\beta_i)}\\ \psi & = & \exp(-C_i(\delta-1)^2-D_i(\tau-1)^2) \end{array}\f$
 ResidualHelmholtzSAFTAssociating   | \f$ \alpha^r = am\left(\ln X-\frac{X}{2}+\frac{1}{2}\right); \f$

--- a/include/MatrixMath.h
+++ b/include/MatrixMath.h
@@ -28,6 +28,43 @@ template<typename T> struct VectorNd<0,T> {


 namespace CoolProp{
+    
+/// Some shortcuts and regularly needed operations
+template<class T> std::size_t         num_rows  (            std::vector<T>   const& in){ return in.size(); }
+template<class T> std::size_t         num_rows  (std::vector<std::vector<T> > const& in){ return in.size(); }
+
+template<class T> std::size_t         max_cols  (std::vector<std::vector<T> > const& in){
+	std::size_t cols = 0;
+	std::size_t col  = 0;
+	for (std::size_t i = 0; i < in.size(); i++) {
+		col = in[i].size();
+		if (cols<col) {cols = col;}
+    }
+	return cols;
+};
+template<class T> bool                is_squared(std::vector<std::vector<T> > const& in){
+	std::size_t cols = max_cols(in);
+	if (cols!=num_rows(in)) { return false;}
+	else {
+		for (std::size_t i = 0; i < in.size(); i++) {
+			if (cols!=in[i].size()) {return false; }
+		}
+	}
+	return true;
+};
+
+template<class T> std::size_t         num_cols  (            std::vector<T>   const& in){ return 1; }
+template<class T> std::size_t         num_cols  (std::vector<std::vector<T> > const& in){
+	if (num_rows(in)>0) {
+		if (is_squared(in)) {
+			return in[0].size();
+		} else {
+			return max_cols(in);
+		}
+	} else {
+		return 0;
+	}
+};    

 /// Convert vectors and matrices
 /** Conversion functions for the different kinds of object-like
@@ -364,42 +401,7 @@ template <class T> std::string mat_to_string(const Eigen::Matrix<T,Eigen::Dynami
 //template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A, const char *fmt);


-/// Some shortcuts and regularly needed operations
-template<class T> std::size_t         num_rows  (            std::vector<T>   const& in){ return in.size(); }
-template<class T> std::size_t         num_rows  (std::vector<std::vector<T> > const& in){ return in.size(); }

-template<class T> std::size_t         max_cols  (std::vector<std::vector<T> > const& in){
-	std::size_t cols = 0;
-	std::size_t col  = 0;
-	for (std::size_t i = 0; i < in.size(); i++) {
-		col = in[i].size();
-		if (cols<col) {cols = col;}
-    }
-	return cols;
-};
-template<class T> bool                is_squared(std::vector<std::vector<T> > const& in){
-	std::size_t cols = max_cols(in);
-	if (cols!=num_rows(in)) { return false;}
-	else {
-		for (std::size_t i = 0; i < in.size(); i++) {
-			if (cols!=in[i].size()) {return false; }
-		}
-	}
-	return true;
-};
-
-template<class T> std::size_t         num_cols  (            std::vector<T>   const& in){ return 1; }
-template<class T> std::size_t         num_cols  (std::vector<std::vector<T> > const& in){
-	if (num_rows(in)>0) {
-		if (is_squared(in)) {
-			return in[0].size();
-		} else {
-			return max_cols(in);
-		}
-	} else {
-		return 0;
-	}
-};


 /*
--- a/src/AbstractState.cpp
+++ b/src/AbstractState.cpp
@@ -338,7 +338,7 @@ TEST_CASE("Check AbstractState","[AbstractState]")
 {
    SECTION("bad backend")
    {
-        CHECK_THROWS(std::tr1::shared_ptr<CoolProp::AbstractState> Water(CoolProp::AbstractState::factory("DEFINITELY_A_BAD_BACKEND", "Water")));
+        CHECK_THROWS(shared_ptr<CoolProp::AbstractState> Water(CoolProp::AbstractState::factory("DEFINITELY_A_BAD_BACKEND", "Water")));
    }
 }

--- a/src/Backends/Helmholtz/ExcessHEFunction.cpp
+++ b/src/Backends/Helmholtz/ExcessHEFunction.cpp
@@ -111,6 +111,7 @@ void ExcessTerm::construct(const std::vector<CoolPropFluid*> &components)

    for (unsigned int i = 0; i < N; ++i)
    {
+        DepartureFunctionMatrix[i].resize(N);
        for (unsigned int j = 0; j < N; ++j)
        {
            if (i == j){ continue; }
@@ -280,7 +281,7 @@ double ExcessTerm::d2alphar_dxi_dDelta(double tau, double delta, const std::vect

 GERG2008DepartureFunction::GERG2008DepartureFunction(const std::vector<double> &n,const std::vector<double> &d,const std::vector<double> &t,
                                                     const std::vector<double> &eta,const std::vector<double> &epsilon,const std::vector<double> &beta,
-                                                     const std::vector<double> &gamma, int Npower)
+                                                     const std::vector<double> &gamma, unsigned int Npower)
 {

    /// Break up into power and gaussian terms
--- a/src/Backends/Helmholtz/ExcessHEFunction.h
+++ b/src/Backends/Helmholtz/ExcessHEFunction.h
@@ -115,7 +115,7 @@ public:
 	GERG2008DepartureFunction(){};
    GERG2008DepartureFunction(const std::vector<double> &n,const std::vector<double> &d,const std::vector<double> &t,
                              const std::vector<double> &eta,const std::vector<double> &epsilon,const std::vector<double> &beta,
-                              const std::vector<double> &gamma, int Npower);
+                              const std::vector<double> &gamma, unsigned int Npower);
 	~GERG2008DepartureFunction(){};
 	double alphar(double tau, double delta);
 	double dalphar_dDelta(double tau, double delta);
--- a/src/Backends/Helmholtz/FlashRoutines.cpp
+++ b/src/Backends/Helmholtz/FlashRoutines.cpp
@@ -69,10 +69,10 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
    }
    else
    {
-        // Set some imput options
+        // Set some input options
        SaturationSolvers::mixture_VLE_IO options;
        options.sstype = SaturationSolvers::imposed_T;
-        options.Nstep_max = 5;
+        options.Nstep_max = 20;

        // Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
        long double pguess = SaturationSolvers::saturation_preconditioner(&HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions);
@@ -81,7 +81,10 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
        pguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess);

        // Actually call the successive substitution solver
-        SaturationSolvers::successive_substitution(&HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options);
+        SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options);
+
+        HEOS._p = options.p;
+        HEOS._rhomolar = 1/(HEOS._Q/options.rhomolar_vap+(1-HEOS._Q)/options.rhomolar_liq);
    }
 }
 void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
@@ -130,10 +133,15 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
        Tguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions, Tguess);

        // Actually call the successive substitution solver
-        SaturationSolvers::successive_substitution(&HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io);
+        SaturationSolvers::successive_substitution(HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io);

-        PhaseEnvelope::PhaseEnvelope_GV ENV_GV;
-        ENV_GV.build(&HEOS, HEOS.mole_fractions, HEOS.K, io);
+        // Load the outputs
+        HEOS._p = HEOS.SatV->p();
+        HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar());
+        HEOS._T = HEOS.SatL->T();
+
+        //PhaseEnvelope::PhaseEnvelope_GV ENV_GV;
+        //ENV_GV.build(&HEOS, HEOS.mole_fractions, HEOS.K, io);
    }
 }
 // D given and one of P,H,S,U
--- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h
+++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h
@@ -311,6 +311,13 @@ protected:
        EOS.reduce.rhomolar = cpjson::get_double(reducing_state,"rhomolar");
        EOS.reduce.p = cpjson::get_double(reducing_state,"p");

+        EOS.sat_min_liquid.T = cpjson::get_double(satminL_state, "T");
+        EOS.sat_min_liquid.p = cpjson::get_double(satminL_state, "p");
+        EOS.sat_min_liquid.rhomolar = cpjson::get_double(satminL_state, "rhomolar");
+        EOS.sat_min_vapor.T = cpjson::get_double(satminV_state, "T");
+        EOS.sat_min_vapor.p = cpjson::get_double(satminV_state, "p");
+        EOS.sat_min_vapor.rhomolar = cpjson::get_double(satminV_state, "rhomolar");
+
        /// \todo: define limits of EOS better
        EOS.limits.Tmin = cpjson::get_double(satminL_state, "T");
        EOS.Ttriple = EOS.limits.Tmin;
--- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
+++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
@@ -50,6 +50,7 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector<CoolPropFluid*> comp

    // Copy the components
    this->components = components;
+    this->N = components.size();

    if (components.size() == 1){
        is_pure_or_pseudopure = true;
@@ -81,14 +82,27 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector<CoolPropFluid*> comp
 }
 void HelmholtzEOSMixtureBackend::set_mole_fractions(const std::vector<long double> &mole_fractions)
 {
-    if (mole_fractions.size() != components.size())
+    if (mole_fractions.size() != N)
    {
-        throw ValueError(format("size of mole fraction vector [%d] does not equal that of component vector [%d]",mole_fractions.size(), components.size()));
+        throw ValueError(format("size of mole fraction vector [%d] does not equal that of component vector [%d]",mole_fractions.size(), N));
+    }
+    // Copy values without reallocating memory
+    this->resize(N);
+    std::copy( mole_fractions.begin(), mole_fractions.end(), this->mole_fractions.begin() );
+    // Resize the vectors for the liquid and vapor,  but only if they are in use
+    if (this->SatL.get() != NULL){
+        this->SatL->resize(N);
+    }
+    if (this->SatV.get() != NULL){
+        this->SatV->resize(N);
    }
-    this->mole_fractions = mole_fractions;
-    this->K.resize(mole_fractions.size());
-    this->lnK.resize(mole_fractions.size());
 };
+void HelmholtzEOSMixtureBackend::resize(unsigned int N)
+{
+    this->mole_fractions.resize(N);
+    this->K.resize(N);
+    this->lnK.resize(N);
+}
 void HelmholtzEOSMixtureBackend::set_reducing_function()
 {
    Reducing.set(ReducingFunction::factory(components));
@@ -1343,7 +1357,9 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double
    try{
        // First we try with Newton's method with analytic derivative
        double rhomolar = Newton(resid, rhomolar_guess, 1e-8, 100, errstring);
-        if (!ValidNumber(rhomolar)){throw ValueError();}
+        if (!ValidNumber(rhomolar)){
+            throw ValueError();
+        }
        return rhomolar;
    }
    catch(std::exception &)
--- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
+++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
@@ -19,14 +19,13 @@ protected:
    std::vector<CoolPropFluid*> components; ///< The components that are in use

    bool is_pure_or_pseudopure; ///< A flag for whether the substance is a pure or pseudo-pure fluid (true) or a mixture (false)
-    std::vector<long double> mole_fractions; ///< The mole fractions of the components
-    std::vector<long double> mole_fractions_liq, ///< The mole fractions of the saturated liquid
-                             mole_fractions_vap, ///< The mole fractions of the saturated vapor
-                             K, ///< The K factors for the components
+    std::vector<long double> mole_fractions; ///< The bulk mole fractions of the mixture
+    std::vector<long double> K, ///< The K factors for the components
                             lnK; ///< The natural logarithms of the K factors of the components

    SimpleState _crit;
    int imposed_phase_index;
+    int N; // Number of components
 public:
    HelmholtzEOSMixtureBackend(){imposed_phase_index = -1;};
    HelmholtzEOSMixtureBackend(std::vector<CoolPropFluid*> components, bool generate_SatL_and_SatV = true);
@@ -46,6 +45,7 @@ public:
    std::vector<long double> &get_K(){return K;};
    std::vector<long double> &get_lnK(){return lnK;};

+    void resize(unsigned int N);
    shared_ptr<HelmholtzEOSMixtureBackend> SatL, SatV; ///<

    void update(long input_pair, double value1, double value2);
@@ -74,7 +74,8 @@ public:
    */
    void set_mole_fractions(const std::vector<long double> &mole_fractions);

-    const std::vector<long double> &get_mole_fractions(){return mole_fractions;};
+    std::vector<long double> &get_mole_fractions(){return mole_fractions;};
+    const std::vector<long double> &get_const_mole_fractions(){return mole_fractions;};

    /// Set the mass fractions
    /**
--- a/src/Backends/Helmholtz/VLERoutines.cpp
+++ b/src/Backends/Helmholtz/VLERoutines.cpp
@@ -10,22 +10,22 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
    /*
    This function is inspired by the method of Akasaka:

-    R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from 
-    Helmholtz Energy Equations of State", 
+    R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
+    Helmholtz Energy Equations of State",
    Journal of Thermal Science and Technology v3 n3,2008

    Ancillary equations are used to get a sensible starting point
    */
    std::vector<long double> negativer(3,_HUGE), v;
    std::vector<std::vector<long double> > J(3, std::vector<long double>(3,_HUGE));
-    
+
    HEOS->calc_reducing_state();
    const SimpleState & reduce = HEOS->get_reducing();
-    shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL, 
+    shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL,
                                           SatV = HEOS->SatV;
    const std::vector<long double> & mole_fractions = HEOS->get_mole_fractions();
    const long double R_u = HEOS->gas_constant();
-    
+
    long double T, rhoL,rhoV;
    long double deltaL=0, deltaV=0, tau=0, error;
    int iter=0;
@@ -35,7 +35,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
    {
        if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_PL || options.specified_variable == saturation_PHSU_pure_options::IMPOSED_PV)
        {
-            // Invert liquid density ancillary to get temperature 
+            // Invert liquid density ancillary to get temperature
            // TODO: fit inverse ancillaries too
            T = HEOS->get_components()[0]->ancillaries.pL.invert(specified_value);
        }
@@ -76,7 +76,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l

        double pL = SatL->p();
        double pV = SatV->p();
-        
+
        // These derivatives are needed for both cases
        double alpharL = SatL->alphar();
        double alpharV = SatV->alphar();
@@ -171,7 +171,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
        v = linsolve(J, negativer);

        tau += options.omega*v[0];
-        
+
        if (options.use_logdelta){
            deltaL = exp(log(deltaL)+options.omega*v[1]);
            deltaV = exp(log(deltaV)+options.omega*v[2]);
@@ -184,7 +184,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
        rhoL = deltaL*reduce.rhomolar;
        rhoV = deltaV*reduce.rhomolar;
        T = reduce.T/tau;
-        
+
        error = sqrt(pow(negativer[0], 2)+pow(negativer[1], 2)+pow(negativer[2], 2));
        iter++;
        if (T < 0)
@@ -202,21 +202,21 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
    /*
    This function is inspired by the method of Akasaka:

-    R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from 
-    Helmholtz Energy Equations of State", 
+    R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
+    Helmholtz Energy Equations of State",
    Journal of Thermal Science and Technology v3 n3,2008

    Ancillary equations are used to get a sensible starting point
    */
    std::vector<long double> r(2,_HUGE), v;
    std::vector<std::vector<long double> > J(2, std::vector<long double>(2,_HUGE));
-    
+
    HEOS->calc_reducing_state();
    const SimpleState & reduce = HEOS->get_reducing();
-    shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL, 
+    shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL,
                                           SatV = HEOS->SatV;
    const std::vector<long double> & mole_fractions = HEOS->get_mole_fractions();
-    
+
    long double T, rhoL,rhoV;
    long double deltaL=0, deltaV=0, tau=0, error;
    int iter=0;
@@ -226,7 +226,7 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
    {
        if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL)
        {
-            // Invert liquid density ancillary to get temperature 
+            // Invert liquid density ancillary to get temperature
            // TODO: fit inverse ancillaries too
            T = HEOS->get_components()[0]->ancillaries.rhoL.invert(rhomolar);
            rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T);
@@ -234,7 +234,7 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
        }
        else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV)
        {
-            // Invert vapor density ancillary to get temperature 
+            // Invert vapor density ancillary to get temperature
            // TODO: fit inverse ancillaries too
            T = HEOS->get_components()[0]->ancillaries.rhoV.invert(rhomolar);
            rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T);
@@ -265,7 +265,7 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long

        double pL = SatL->p();
        double pV = SatV->p();
-        
+
        // These derivatives are needed for both cases
        double dalphar_dtauL = SatL->dalphar_dTau();
        double dalphar_dtauV = SatV->dalphar_dTau();
@@ -275,8 +275,8 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
        double alpharV = SatV->alphar();
        double dalphar_ddeltaL = SatL->dalphar_dDelta();
        double dalphar_ddeltaV = SatV->dalphar_dDelta();
-        
-        
+
+
        // -r_1
        r[0] = -(deltaV*(1+deltaV*dalphar_ddeltaV)-deltaL*(1+deltaL*dalphar_ddeltaL));
        // -r_2
@@ -340,7 +340,7 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
        rhoL = deltaL*reduce.rhomolar;
        rhoV = deltaV*reduce.rhomolar;
        T = reduce.T/tau;
-        
+
        error = sqrt(pow(r[0], 2)+pow(r[1], 2));
        iter++;
        if (T < 0)
@@ -367,19 +367,19 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
    // Start with the method of Akasaka

    /*
-    This function implements the method of Akasaka 
+    This function implements the method of Akasaka

-    R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from 
-    Helmholtz Energy Equations of State", 
+    R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
+    Helmholtz Energy Equations of State",
    Journal of Thermal Science and Technology v3 n3,2008

    Ancillary equations are used to get a sensible starting point
    */
-    
+
    HEOS->calc_reducing_state();
    const SimpleState & reduce = HEOS->get_reducing();
    long double R_u = HEOS->calc_gas_constant();
-    shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL, 
+    shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL,
                                           SatV = HEOS->SatV;

    long double rhoL,rhoV,JL,JV,KL,KV,dJL,dJV,dKL,dKV;
@@ -390,7 +390,7 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
    {
        if (options.use_guesses)
        {
-            // Use the guesses provided in the options structure            
+            // Use the guesses provided in the options structure
            rhoL = options.rhoL;
            rhoV = options.rhoV;
        }
@@ -405,7 +405,7 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
                rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T);
                rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T);
            }
-        }    
+        }

        // Apply a single step of Newton's method to improve guess value for liquid
        // based on the error between the gas pressure (which is usually very close already)
@@ -451,7 +451,7 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
        double dalphar_ddeltaV = SatV->dalphar_dDelta();
        double d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
        double d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
-        
+
        JL = deltaL * (1 + deltaL*dalphar_ddeltaL);
        JV = deltaV * (1 + deltaV*dalphar_ddeltaV);
        KL = deltaL*dalphar_ddeltaL + alpharL + log(deltaL);
@@ -459,13 +459,13 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE

        PL = R_u*reduce.rhomolar*T*JL;
        PV = R_u*reduce.rhomolar*T*JV;
-        
+
        // At low pressure, the magnitude of d2alphar_ddelta2L and d2alphar_ddelta2V are enormous, truncation problems arise for all the partials
        dJL = 1 + 2*deltaL*dalphar_ddeltaL + deltaL*deltaL*d2alphar_ddelta2L;
        dJV = 1 + 2*deltaV*dalphar_ddeltaV + deltaV*deltaV*d2alphar_ddelta2V;
        dKL = 2*dalphar_ddeltaL + deltaL*d2alphar_ddelta2L + 1/deltaL;
        dKV = 2*dalphar_ddeltaV + deltaV*d2alphar_ddelta2V + 1/deltaV;
-        
+
        DELTA = dJV*dKL-dJL*dKV;

        error = fabs(KL-KV)+fabs(JL-JV);
@@ -473,9 +473,9 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
        //  Get the predicted step
        stepL = options.omega/DELTA*( (KV-KL)*dJV-(JV-JL)*dKV);
        stepV = options.omega/DELTA*( (KV-KL)*dJL-(JV-JL)*dKL);
-        
+
        if (deltaL+stepL > 1 && deltaV+stepV < 1 && deltaV+stepV > 0){
-            deltaL += stepL;  deltaV += stepV; 
+            deltaL += stepL;  deltaV += stepV;
            rhoL = deltaL*reduce.rhomolar; rhoV = deltaV*reduce.rhomolar;
        }
        else{
@@ -497,62 +497,67 @@ void SaturationSolvers::x_and_y_from_K(long double beta, const std::vector<long
        x[i] = z[i]/denominator;
        y[i] = K[i]*z[i]/denominator;
    }
-    //normalize_vector(x);
-    //normalize_vector(y);
 }

-long double SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend *HEOS, const long double beta, long double T, long double p, const std::vector<long double> &z, 
+long double SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS, const long double beta, long double T, long double p, const std::vector<long double> &z,
                                                       std::vector<long double> &K, mixture_VLE_IO &options)
 {
    int iter = 1;
    long double change, f, df, deriv_liq, deriv_vap;
    std::size_t N = z.size();
-    std::vector<long double> x, y, ln_phi_liq, ln_phi_vap;
-    ln_phi_liq.resize(N); ln_phi_vap.resize(N); x.resize(N); y.resize(N);
+    std::vector<long double> ln_phi_liq, ln_phi_vap;
+    ln_phi_liq.resize(N); ln_phi_vap.resize(N);

+    std::vector<long double> &x = HEOS.SatL->get_mole_fractions(), &y = HEOS.SatV->get_mole_fractions();
    x_and_y_from_K(beta, K, z, x, y);
-    shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())), 
-                                           SatV(new HelmholtzEOSMixtureBackend(HEOS->get_components()));
-    SatL->specify_phase(iphase_liquid);
-    SatV->specify_phase(iphase_gas);

-    SatL->set_mole_fractions(x);
-    SatV->set_mole_fractions(y);
-    long double rhomolar_liq = SatL->solver_rho_Tp_SRK(T, p, iphase_liquid); // [mol/m^3]
-    long double rhomolar_vap = SatV->solver_rho_Tp_SRK(T, p, iphase_gas); // [mol/m^3]
-    
+    HEOS.SatL->specify_phase(iphase_liquid);
+    HEOS.SatV->specify_phase(iphase_gas);
+
+    normalize_vector(x);
+    normalize_vector(y);
+
+    HEOS.SatL->set_mole_fractions(x);
+    HEOS.SatV->set_mole_fractions(y);
+    HEOS.SatL->calc_reducing_state();
+    HEOS.SatV->calc_reducing_state();
+    long double rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(T, p, iphase_liquid); // [mol/m^3]
+    long double rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(T, p, iphase_gas); // [mol/m^3]
+
+    HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq*1.3);
+    HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap);
+
    do
    {
-        
-        SatL->update_TP_guessrho(T, p, rhomolar_liq);
-        SatV->update_TP_guessrho(T, p, rhomolar_vap);
+        HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
+        HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());

        f = 0;
        df = 0;

        for (std::size_t i = 0; i < N; ++i)
        {
-            ln_phi_liq[i] = SatL->mixderiv_ln_fugacity_coefficient(i);
-            ln_phi_vap[i] = SatV->mixderiv_ln_fugacity_coefficient(i);
+            ln_phi_liq[i] = HEOS.SatL->mixderiv_ln_fugacity_coefficient(i);
+            ln_phi_vap[i] = HEOS.SatV->mixderiv_ln_fugacity_coefficient(i);

            if (options.sstype == imposed_p){
-                deriv_liq = SatL->mixderiv_dln_fugacity_coefficient_dT__constp_n(i);
-                deriv_vap = SatV->mixderiv_dln_fugacity_coefficient_dT__constp_n(i);
+                deriv_liq = HEOS.SatL->mixderiv_dln_fugacity_coefficient_dT__constp_n(i);
+                deriv_vap = HEOS.SatV->mixderiv_dln_fugacity_coefficient_dT__constp_n(i);
            }
            else if (options.sstype == imposed_T){
-                deriv_liq = SatL->mixderiv_dln_fugacity_coefficient_dp__constT_n(i);
-                deriv_vap = SatV->mixderiv_dln_fugacity_coefficient_dp__constT_n(i);
+                deriv_liq = HEOS.SatL->mixderiv_dln_fugacity_coefficient_dp__constT_n(i);
+                deriv_vap = HEOS.SatV->mixderiv_dln_fugacity_coefficient_dp__constT_n(i);
            }
            else {throw ValueError();}

            K[i] = exp(ln_phi_liq[i]-ln_phi_vap[i]);
-            
+
            f += z[i]*(K[i]-1)/(1-beta+beta*K[i]);

            double dfdK = K[i]*z[i]/pow(1-beta+beta*K[i],(int)2);
            df += dfdK*(deriv_liq-deriv_vap);
        }
-        
+
        change = -f/df;

        if (options.sstype == imposed_p){
@@ -563,8 +568,8 @@ long double SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBacken
        }

        x_and_y_from_K(beta, K, z, x, y);
-        SatL->set_mole_fractions(x);
-        SatV->set_mole_fractions(y);
+        HEOS.SatL->set_mole_fractions(x);
+        HEOS.SatV->set_mole_fractions(y);

        iter += 1;
        if (iter > 50)
@@ -572,32 +577,24 @@ long double SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBacken
            return _HUGE;
            //throw ValueError(format("saturation_p was unable to reach a solution within 50 iterations"));
        }
-        rhomolar_liq = SatL->rhomolar();
-        rhomolar_vap = SatV->rhomolar();
    }
    while(fabs(f) > 1e-12 || iter < options.Nstep_max);

-    SatL->update_TP_guessrho(T, p, rhomolar_liq);
-    SatV->update_TP_guessrho(T, p, rhomolar_vap);
-    
-    double pL = SatL->calc_pressure();
-    double pV = SatV->calc_pressure();
+    HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq);
+    HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap);

-
-    options.rhomolar_liq = SatL->rhomolar();
-    options.rhomolar_vap = SatV->rhomolar();
-    options.p = pL;
-    options.T = T;
-    options.x = SatL->get_mole_fractions();
-    options.y = SatV->get_mole_fractions();
+    options.p = HEOS.SatL->p();
+    options.T = HEOS.SatL->T();
+    options.rhomolar_liq = HEOS.SatL->rhomolar();
+    options.rhomolar_vap = HEOS.SatV->rhomolar();
 }

 void SaturationSolvers::newton_raphson_VLE_GV::resize(unsigned int N)
 {
    this->N = N;
-    x.resize(N); 
-    y.resize(N); 
-    phi_ij_liq.resize(N); 
+    x.resize(N);
+    y.resize(N);
+    phi_ij_liq.resize(N);
    phi_ij_vap.resize(N);
    dlnphi_drho_liq.resize(N);
    dlnphi_drho_vap.resize(N);
@@ -621,7 +618,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtur
    std::size_t N = K.size();
    resize(N);

-    shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())), 
+    shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())),
                                           SatV(new HelmholtzEOSMixtureBackend(HEOS->get_components()));
    SatL->specify_phase(iphase_liquid);
    SatV->specify_phase(iphase_gas);
@@ -633,7 +630,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtur

    // Build the Jacobian and residual vectors for given inputs of K_i,T,p
    build_arrays(HEOS,beta0,T0,rhomolar_liq0,rhomolar_vap0,z,K);
-    
+
    // Make copies of the base
    std::vector<long double> r0 = r;
    STLMatrix J0 = J;
@@ -644,7 +641,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtur
        for (std::size_t j = 0; j < N+2; ++j)
        {
            Jnum[i][j] = _HUGE;
-        }    
+        }
    }

    for (std::size_t j = 0; j < N; ++j)
@@ -660,7 +657,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtur
        std::cout << vec_to_string(get_col(Jnum,j),"%12.11f") << std::endl;
        std::cout << vec_to_string(get_col(J,j),"%12.11f") << std::endl;
    }
-    
+
    build_arrays(HEOS,beta0,T0+1e-6,rhomolar_liq0,rhomolar_vap0,z,K);
    std::vector<long double> r1 = r, JN = r;
    for (std::size_t i = 0; i < r.size(); ++i)
@@ -689,9 +686,9 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend *
    pre_call();
    resize(K.size());

-    shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())), 
+    shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())),
                                           SatV(new HelmholtzEOSMixtureBackend(HEOS->get_components()));
-    SatL->specify_phase(iphase_liquid); // So it will always just use single-phase solution 
+    SatL->specify_phase(iphase_liquid); // So it will always just use single-phase solution
    SatV->specify_phase(iphase_gas); // So it will always just use single-phase solution

    do
@@ -699,7 +696,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend *
        // Build the Jacobian and residual vectors for given inputs of K_i,T,p
        build_arrays(HEOS,IO.beta,IO.T,IO.rhomolar_liq,IO.rhomolar_vap,z,K);

-        // Solve for the step; v is the step with the contents 
+        // Solve for the step; v is the step with the contents
        // [delta(lnK0), delta(lnK1), ..., delta(lnT), delta(lnrho')]
        std::vector<long double> v = linsolve(J, negative_r);

@@ -723,15 +720,15 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend *
            std::cout << "nr = " << vec_to_string(r,"%16.15g");*/
            throw ValueError("Temperature or p has bad value");
        }
-        
+
        //std::cout << iter << " " << T << " " << p << " " << error_rms << std::endl;
        iter++;
    }
    while(this->error_rms > 1e-8 && max_rel_change > 1000*LDBL_EPSILON && iter < IO.Nstep_max);
    Nsteps = iter;
    IO.p = p;
-    IO.x = x; // Mole fractions in liquid
-    IO.y = y; // Mole fractions in vapor
+    IO.x = &x; // Mole fractions in liquid
+    IO.y = &y; // Mole fractions in vapor
 }

 void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureBackend *HEOS, long double beta, long double T, long double rhomolar_liq, const long double rhomolar_vap, const std::vector<long double> &z, std::vector<long double> &K)
@@ -778,11 +775,11 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
            // I think this is wrong.
            phi_ij_vap[j] = SatV->mixderiv_ndln_fugacity_coefficient_dnj__constT_p(i,j) + (SatV->mixderiv_partial_molar_volume(i)/(SatV->gas_constant()*T)-1/p)*SatV->mixderiv_ndpdni__constT_V_nj(i); ; // 7.126 from GERG monograph
        }
-        
+
        r[i] = log(K[i]) + ln_phi_vap - ln_phi_liq;
        // dF_i/d(ln(K_j))
        for (unsigned int j = 0; j < N; ++j)
-        {	
+        {
            J[i][j] = K[j]*z[j]/pow(1-beta+beta*K[j],(int)2)*((1-beta)*phi_ij_vap[j]+beta*phi_ij_liq[j])+Kronecker_delta(i,j);
        }
        // dF_{i}/d(ln(T))
@@ -794,10 +791,10 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
    double summer1 = 0;
    for (unsigned int i = 0; i < N; ++i)
    {
-        // Although the definition of this term is given by 
-        // y[i]-x[i], when x and y are normalized, you get 
+        // Although the definition of this term is given by
+        // y[i]-x[i], when x and y are normalized, you get
        // the wrong values.  Why? No idea.
-        summer1 += z[i]*(K[i]-1)/(1-beta+beta*K[i]); 
+        summer1 += z[i]*(K[i]-1)/(1-beta+beta*K[i]);
    }
    r[N] = summer1;

@@ -811,7 +808,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
    // For the residual term F_{N+1} = p'-p''
    r[N+1] = p_liq-p_vap;
    for (unsigned int j = 0; j < N; ++j)
-    {	
+    {
        J[N+1][j] = HEOS->gas_constant()*T*K[j]*z[j]/pow(1-beta+beta*K[j],(int)2)*((1-beta)*dlnphi_drho_vap[j]+beta*dlnphi_drho_liq[j]);
    }
    // dF_{N+1}/d(ln(T))
@@ -837,7 +834,7 @@ void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend *HEOS, co
    SaturationSolvers::mixture_VLE_IO IO_NRVLE = IO;
    bubble.resize(z.size());
    dew.resize(z.size());
-        
+
    // HACK
    IO_NRVLE.beta = 1.0;
    IO_NRVLE.Nstep_max = 30;
@@ -881,7 +878,7 @@ void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend *HEOS, co
            NRVLE.check_Jacobian(HEOS,z,K,IO_NRVLE);
        }*/
        NRVLE.call(HEOS,z,K,IO_NRVLE);
-        dew.store_variables(IO_NRVLE.T,IO_NRVLE.p,IO_NRVLE.rhomolar_liq,IO_NRVLE.rhomolar_vap,K,IO_NRVLE.x,IO_NRVLE.y);
+        dew.store_variables(IO_NRVLE.T,IO_NRVLE.p,IO_NRVLE.rhomolar_liq,IO_NRVLE.rhomolar_vap,K,*(IO_NRVLE.x),*(IO_NRVLE.y));
        iter ++;
        std::cout << format("%g %g %g %g %g %d %g\n",IO_NRVLE.p,IO_NRVLE.rhomolar_liq,IO_NRVLE.rhomolar_vap,IO_NRVLE.T,K[0],NRVLE.Nsteps,factor);
        if (iter < 5){continue;}
@@ -902,4 +899,4 @@ void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend *HEOS, co
    }
 }

-} /* namespace CoolProp*/
+} /* namespace CoolProp*/
--- a/src/Backends/Helmholtz/VLERoutines.h
+++ b/src/Backends/Helmholtz/VLERoutines.h
@@ -31,7 +31,7 @@ namespace SaturationSolvers
    {
        int sstype, Nstep_max;
        long double rhomolar_liq, rhomolar_vap, p, T, beta;
-        std::vector<long double> x,y,K;
+        std::vector<long double> *x, *y, *K;
    };

    /*! Returns the natural logarithm of K for component i using the method from Wilson as in
@@ -67,7 +67,13 @@ namespace SaturationSolvers
    */
    void saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, long double specified_value, saturation_PHSU_pure_options &options);

-    long double successive_substitution(HelmholtzEOSMixtureBackend *HEOS, const long double beta, long double T, long double p, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &options);
+    long double successive_substitution(HelmholtzEOSMixtureBackend &HEOS,
+                                        const long double beta,
+                                        long double T,
+                                        long double p,
+                                        const std::vector<long double> &z,
+                                        std::vector<long double> &K,
+                                        mixture_VLE_IO &options);
    void x_and_y_from_K(long double beta, const std::vector<long double> &K, const std::vector<long double> &z, std::vector<long double> &x, std::vector<long double> &y);

    /*! A wrapper function around the residual to find the initial guess for the bubble point temperature
@@ -116,9 +122,9 @@ namespace SaturationSolvers
 	    {
            EquationOfState *EOS = (HEOS->get_components())[i]->pEOS; 

-		    ptriple += EOS->ptriple*z[i];
+		    ptriple += EOS->sat_min_liquid.p*z[i];
            pcrit += EOS->reduce.p*z[i];
-		    Ttriple += EOS->Ttriple*z[i];
+		    Ttriple += EOS->sat_min_liquid.T*z[i];
 		    Tcrit += EOS->reduce.T*z[i];
 	    }

--- a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp
+++ b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp
@@ -1233,12 +1233,12 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
            // Skip fluids not in REFPROP
            if (RPName.find("N/A") == 0){continue;}

-            std::tr1::shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
+            ::shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
            double Tr = S1->T_critical();
            S1->update(CoolProp::QT_INPUTS, 0, Tr*0.9);
            double rho_CP = S1->rhomolar();

-            std::tr1::shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
+            ::shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
            S2->update(CoolProp::QT_INPUTS, 0, Tr*0.9);
            double rho_RP = S2->rhomolar();

@@ -1263,12 +1263,12 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
            // Skip fluids not in REFPROP
            if (RPName.find("N/A") == 0){continue;}

-            std::tr1::shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
+            shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
            double Tr = S1->T_critical();
            S1->update(CoolProp::QT_INPUTS, 0, Tr*0.9);
            double cp_CP = S1->cpmolar();

-            std::tr1::shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
+            shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
            S2->update(CoolProp::QT_INPUTS, 0, Tr*0.9);
            double cp_RP = S2->cpmolar();

@@ -1294,13 +1294,13 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
            // Skip fluids not in REFPROP
            if (RPName.find("N/A") == 0){continue;}

-            std::tr1::shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
+            shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
            double Tr = S1->T_critical();
            S1->update(CoolProp::QT_INPUTS, 0, 0.9*Tr);
            double h_CP = S1->hmass();
            double s_CP = S1->smass();

-            std::tr1::shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
+            shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
            S2->update(CoolProp::QT_INPUTS, 0, 0.9*Tr);
            double h_RP = S2->hmass();
            double s_RP = S2->smass();
--- a/src/Backends/REFPROP/REFPROPMixtureBackend.h
+++ b/src/Backends/REFPROP/REFPROPMixtureBackend.h
@@ -49,6 +49,13 @@ public:
 				double value2
 				);

+    long double calc_molar_mass(void){
+        if (!_molar_mass){
+            calc_molar_mass();
+        }
+        return _molar_mass.pt();
+    };
+
 	/// Returns true if REFPROP is supported on this platform
 	bool REFPROP_supported(void);

--- a/src/Helmholtz.cpp
+++ b/src/Helmholtz.cpp
@@ -1722,7 +1722,7 @@ class HelmholtzConsistencyFixture
 public:
    long double numerical, analytic;

-    std::tr1::shared_ptr<CoolProp::BaseHelmholtzTerm> Lead, LogTau, IGPower, PlanckEinstein,
+    shared_ptr<CoolProp::BaseHelmholtzTerm> Lead, LogTau, IGPower, PlanckEinstein,
        CP0Constant, CP0PolyT, Gaussian, Lemmon2005, Power, SAFT, NonAnalytic, Exponential, GERG2008;

    HelmholtzConsistencyFixture(){
@@ -1942,7 +1942,7 @@ std::string derivs[] = {"dTau","dTau2","dTau3","dDelta","dDelta2","dDelta3","dDe

 TEST_CASE_METHOD(HelmholtzConsistencyFixture, "Helmholtz energy derivatives", "[helmholtz]")
 {
-    std::tr1::shared_ptr<CoolProp::BaseHelmholtzTerm> term;
+    shared_ptr<CoolProp::BaseHelmholtzTerm> term;
    std::size_t n = sizeof(terms)/sizeof(terms[0]);
    for (std::size_t i = 0; i < n; ++i)
    {
--- a/src/MatrixMath.cpp
+++ b/src/MatrixMath.cpp
@@ -12,389 +12,181 @@

 namespace CoolProp{

-///*
-//Owe a debt of gratitude to http://sole.ooz.ie/en - very clear treatment of GJ
-//*/
-//template<typename T> void swap_rows(std::vector<std::vector<T> > *A, size_t row1, size_t row2)
-//{
-//	for (size_t col = 0; col < (*A)[0].size(); col++){
-//		std::swap((*A)[row1][col],(*A)[row2][col]);
-//	}
-//}
-//template<typename T> void subtract_row_multiple(std::vector<std::vector<T> > *A, size_t row, T multiple, size_t pivot_row)
-//{
-//	for (size_t col = 0; col < (*A)[0].size(); col++){
-//		(*A)[row][col] -= multiple*(*A)[pivot_row][col];
-//	}
-//}
-//template<typename T> void divide_row_by(std::vector<std::vector<T> > *A, size_t row, T value)
-//{
-//	for (size_t col = 0; col < (*A)[0].size(); col++){
-//		(*A)[row][col] /= value;
-//	}
-//}
-//
-//template<typename T> size_t get_pivot_row(std::vector<std::vector<T> > *A, size_t col)
-//{
-//	int index = col;
-//	T max = 0, val;
-//
-//	for (size_t row = col; row < (*A).size(); row++)
-//	{
-//		val = (*A)[row][col];
-//		if (fabs(val) > max)
-//		{
-//			max = fabs(val);
-//			index = row;
-//		}
-//	}
-//	return index;
-//}
-//
-//
-//template<typename T> std::vector<std::vector<T> > linsolve_Gauss_Jordan(std::vector<std::vector<T> > const& A, std::vector<std::vector<T> > const& B) {
-//	std::vector<std::vector<T> > AB;
-//	std::vector<std::vector<T> > X;
-//	size_t pivot_row;
-//	T pivot_element;
-//
-//	size_t NrowA = num_rows(A);
-//	size_t NrowB = num_rows(B);
-//	size_t NcolA = num_cols(A);
-//	size_t NcolB = num_cols(B);
-//
-//	if (NrowA!=NrowB) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ",NrowA,NrowB));
-//
-//	AB.resize(NrowA, std::vector<T>(NcolA+NcolB, 0));
-//	 X.resize(NrowA, std::vector<T>(NcolB, 0));
-//
-//	// Build the augmented matrix
-//	for (size_t row = 0; row < NrowA; row++){
-//		for (size_t col  = 0; col < NcolA; col++){
-//			AB[row][col] = A[row][col];
-//		}
-//		for (size_t col  = NcolA; col < NcolA+NcolB; col++){
-//			AB[row][col] = B[row][col-NcolA];
-//		}
-//	}
-//
-//	for (size_t col = 0; col < NcolA; col++){
-//		// Find the pivot value
-//		pivot_row = get_pivot_row(&AB, col);
-//
-//		if (fabs(AB[pivot_row][col]) < 10*DBL_EPSILON){ throw ValueError(format("Zero occurred in row %d, the matrix is singular. ",pivot_row));}
-//
-//		if (pivot_row>=col){
-//			// Swap pivot row and current row
-//			swap_rows(&AB, col, pivot_row);
-//		}
-//		// Get the pivot element
-//		pivot_element = AB[col][col];
-//		// Divide the pivot row by the pivot element
-//		divide_row_by(&AB,col,pivot_element);
-//
-//		if (col < NrowA-1)
-//		{
-//			// All the rest of the rows, subtract the value of the [r][c] combination
-//			for (size_t row = col + 1; row < NrowA; row++)
-//			{
-//				subtract_row_multiple(&AB,row,AB[row][col],col);
-//			}
-//		}
-//	}
-//	for (int col = NcolA - 1; col > 0; col--)
-//	{
-//		for (int row = col - 1; row >=0; row--)
-//		{
-//			subtract_row_multiple(&AB,row,AB[row][col],col);
-//		}
-//	}
-//	// Set the output value
-//	for (size_t row = 0; row < NrowA; row++){
-//		for (size_t col  = 0; col < NcolB; col++){
-//			X[row][col] = AB[row][NcolA+col];
-//		}
-//	}
-//	return X;
-//}
-//
-//
-////std::vector<std::vector<double> > linsolve_Gauss_Jordan_reimpl(std::vector<std::vector<double> > const& A, std::vector<std::vector<double> > const& B) {
-////	std::vector<std::vector<double> > AB;
-////	std::vector<std::vector<double> > X;
-////	size_t pivot_row;
-////	double pivot_element;
-////	double tmp_element;
-////
-////	size_t NrowA = num_rows(A);
-////	size_t NrowB = num_rows(B);
-////	size_t NcolA = num_cols(A);
-////	size_t NcolB = num_cols(B);
-////
-////	if (NrowA!=NrowB) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ",NrowA,NrowB));
-////
-////	AB.resize(NrowA, std::vector<double>(NcolA+NcolB, 0));
-////	 X.resize(NrowA, std::vector<double>(NcolB, 0));
-////
-////	// Build the augmented matrix
-////	for (size_t row = 0; row < NrowA; row++){
-////		for (size_t col  = 0; col < NcolA; col++){
-////			AB[row][col] = A[row][col];
-////		}
-////		for (size_t col  = NcolA; col < NcolA+NcolB; col++){
-////			AB[row][col] = B[row][col-NcolA];
-////		}
-////	}
-////
-////	for (size_t col = 0; col < NcolA; col++){
-////		// Find the pivot row
-////		pivot_row     = 0;
-////		pivot_element = 0.0;
-////		for (size_t row = col; row < NrowA; row++){
-////			tmp_element = fabs(AB[row][col]);
-////			if (tmp_element>pivot_element) {
-////				pivot_element = tmp_element;
-////				pivot_row = row;
-////			}
-////		}
-////		// Check for errors
-////		if (AB[pivot_row][col]<1./_HUGE) throw ValueError(format("Zero occurred in row %d, the matrix is singular. ",pivot_row));
-////		// Swap the rows
-////		if (pivot_row>col) {
-////			for (size_t colInt = 0; colInt < NcolA; colInt++){
-////				std::swap(AB[pivot_row][colInt],AB[pivot_row][colInt]);
-////			}
-////		}
-////		// Process the entries below current element
-////		for (size_t row = col; row < NrowA; row++){
-////			// Entries to the right of current element (until end of A)
-////			for (size_t colInt = col+1; colInt < NcolA; colInt++){
-////				// All entries in augmented matrix
-////				for (size_t colFull = col; colFull < NcolA+NcolB; colFull++){
-////					AB[colInt][colFull] -= AB[col][colFull] * AB[colInt][col] / AB[col][col];
-////				}
-////				AB[colInt][col] = 0.0;
-////			}
-////		}
-////	}
-////	return AB;
-////}
-//
-//
-//
-//
-//
-//
-//template<class T> std::vector<std::vector<T> > linsolve(std::vector<std::vector<T> > const& A, std::vector<std::vector<T> > const& B){
-//	return linsolve_Gauss_Jordan(A, B);
-//}
-//
-//template<class T> std::vector<T>               linsolve(std::vector<std::vector<T> > const& A,             std::vector<T>   const& b){
-//	std::vector<std::vector<T> > B;
-//	for (size_t i = 0; i < b.size(); i++){
-//		B.push_back(std::vector<T>(1,b[i]));
-//	}
-//	B = linsolve(A, B);
-//	B[0].resize(B.size(),0.0);
-//	for (size_t i = 1; i < B.size(); i++){
-//		B[0][i] = B[i][0];
-//	}
-//	return B[0];
-//}
-//
-//
-///// Some shortcuts and regularly needed operations
-//template<class T> std::size_t         num_rows  (std::vector<std::vector<T> > const& in){ return in.size(); }
-//template<class T> std::size_t         num_cols  (std::vector<std::vector<T> > const& in){
-//	if (num_rows(in)>0) {
-//		if (is_squared(in)) {
-//			return in[0].size();
-//		} else {
-//			return max_cols(in);
-//		}
-//	} else {
-//		return 0;
-//	}
-//}
-//template<class T> std::size_t         max_cols  (std::vector<std::vector<T> > const& in){
-//	std::size_t cols = 0;
-//	std::size_t col  = 0;
-//	for (std::size_t i = 0; i < in.size(); i++) {
-//		col = in[i].size();
-//		if (cols<col) {cols = col;}
-//    }
-//	return cols;
-//}
-//template<class T> std::vector<T> get_row(std::vector< std::vector<T> > const& in, size_t row) { return in[row]; }
-//template<class T> std::vector<T> get_col(std::vector< std::vector<T> > const& in, size_t col) {
-//	std::size_t sizeX  = in.size();
-//	if (sizeX<1) throw ValueError(format("You have to provide values, a vector length of %d is not valid. ",sizeX));
-//	size_t sizeY = in[0].size();
-//	if (sizeY<1) throw ValueError(format("You have to provide values, a vector length of %d is not valid. ",sizeY));
-//	std::vector<T> out;
-//	for (std::size_t i = 0; i < sizeX; i++) {
-//		sizeY = in[i].size();
-//		if (sizeY-1<col) throw ValueError(format("Your matrix does not have enough entries in row %d, last index %d is less than %d. ",i,sizeY-1,col));
-//		out.push_back(in[i][col]);
-//	}
-//	return out;
-//}
-//template<class T> bool                is_squared(std::vector<std::vector<T> > const& in){
-//	std::size_t cols = max_cols(in);
-//	if (cols!=num_rows(in)) { return false;}
-//	else {
-//		for (std::size_t i = 0; i < in.size(); i++) {
-//			if (cols!=in[i].size()) {return false; }
-//		}
-//	}
-//	return true;
-//}
-//template<class T> std::vector<std::vector<T> > make_squared(std::vector<std::vector<T> > const& in){
-//	std::size_t cols   = max_cols(in);
-//	std::size_t rows   = num_rows(in);
-//	std::size_t maxVal = 0;
-//	std::vector<std::vector<T> > out;
-//	            std::vector<T>   tmp;
-//
-//	if (cols>rows) {maxVal = cols; }
-//	else           {maxVal = rows; }
-//	out.clear();
-//	for (std::size_t i = 0; i < in.size(); i++) {
-//		tmp.clear();
-//		for (std::size_t j = 0; j < in[i].size(); j++) {
-//			tmp.push_back(in[i][j]);
-//		}
-//		while (maxVal>tmp.size()) {
-//			tmp.push_back(0.0);
-//		}
-//		out.push_back(tmp);
-//    }
-//	// Check rows
-//	tmp.clear();
-//	tmp.resize(maxVal,0.0);
-//	while (maxVal>out.size()) {
-//		out.push_back(tmp);
-//	}
-//	return out;
-//}
-//
-//template<class T>                         T    multiply(            std::vector<T>   const& a,             std::vector<T>   const& b){
-//    return dot_product(a,b);
-//
-//}
-//template<class T>             std::vector<T>   multiply(std::vector<std::vector<T> > const& A,             std::vector<T>   const& b){
-//	std::vector<std::vector<T> > B;
-//	for (size_t i = 0; i < b.size(); i++){
-//		B.push_back(std::vector<T>(1,b[i]));
-//	}
-//	B = multiply(A, B);
-//	B[0].resize(B.size(),0.0);
-//	for (size_t i = 1; i < B.size(); i++){
-//		B[0][i] = B[i][0];
-//	}
-//	return B[0];
-//}
-//
-//template<class T> std::vector<std::vector<T> > multiply(std::vector<std::vector<T> > const& A, std::vector<std::vector<T> > const& B){
-//	if (num_cols(A) != num_rows(B)){
-//		throw ValueError(format("You have to provide matrices with the same columns and rows: %d is not equal to %d. ",num_cols(A),num_rows(B)));
-//	}
-//	size_t rows = num_rows(A);
-//	size_t cols = num_cols(B);
-//	T tmp;
-//	std::vector<std::vector<T> > outVec;
-//		        std::vector<T>   tmpVec;
-//    outVec.clear();
-//	for (size_t i = 0; i < rows; i++){
-//		tmpVec.clear();
-//		for (size_t j = 0; j < cols; j++){
-//			tmp = 0.0;
-//			for (size_t k = 0; k < num_cols(A); k++){
-//				tmp += A[i][k] * B[k][j];
-//			}
-//			tmpVec.push_back(tmp);
-//		}
-//		outVec.push_back(tmpVec);
-//	}
-//	return outVec;
-//}
-//
-//template<class T> T              dot_product(std::vector<T> const& a, std::vector<T> const& b){
-//	if (a.size()==b.size()){
-//		return std::inner_product(a.begin(), a.end(), b.begin(), 0.0);
-//	}
-//	throw ValueError(format("You have to provide vectors with the same length: %d is not equal to %d. ",a.size(),b.size()));
-//}
-//
-//template<class T> std::vector<T> cross_product(std::vector<T> const& a, std::vector<T> const& b){
-//	throw NotImplementedError("The cross product function has not been implemented, yet");
-//}
-//
-//template<class T> std::vector< std::vector<T> > transpose(std::vector<std::vector<T> > const& in){
-//	size_t sizeX = in.size();
-//	if (sizeX<1) throw ValueError(format("You have to provide values, a vector length of %d is not a valid. ",sizeX));
-//	size_t sizeY    = in[0].size();
-//	size_t sizeYOld = sizeY;
-//	if (sizeY<1) throw ValueError(format("You have to provide values, a vector length of %d is not a valid. ",sizeY));
-//	std::vector< std::vector<T> > out(sizeY,std::vector<T>(sizeX));
-//	for (size_t i = 0; i < sizeX; ++i){
-//		sizeY = in[i].size();
-//		if (sizeY!=sizeYOld) throw ValueError(format("You have to provide a rectangular matrix: %d is not equal to %d. ",sizeY,sizeYOld));
-//		for (size_t j = 0; j < sizeY; ++j){
-//			out[j][i] = in[i][j];
-//		}
-//	}
-//	return out;
-//}
-//
-//template<class T> std::vector< std::vector<T> >    invert(std::vector<std::vector<T> > const& in){
-//	if (!is_squared(in)) throw ValueError(format("Only square matrices can be inverted: %d is not equal to %d. ",num_rows(in),num_cols(in)));
-//	std::vector<std::vector<T> > identity;
-//	// Build the identity matrix
-//	size_t dim = num_rows(in);
-//	identity.resize(dim, std::vector<T>(dim, 0));
-//	for (size_t row = 0; row < dim; row++){
-//		identity[row][row] = 1.0;
-//	}
-//	return linsolve(in,identity);
-//}
-//
-//template<class T> std::string vec_to_string(                        T    const& a){
-//	std::stringstream out;
-//	out << format("[ %7.3f ]",a);
-//	return out.str();
-//}
-//
-//template<class T> std::string vec_to_string(            std::vector<T>   const& a) {
-//	return vec_to_string(a,"%7.3g");
-//}
-//template<class T> std::string vec_to_string(            std::vector<T>   const& a, const char *fmt) {
-//	if (a.size()<1) {
-//		return std::string("");
-//	} else {
-//		std::stringstream out;
-//		out << format("[ ");
-//		out << format(fmt,a[0]);
-//		for (size_t j = 1; j < a.size(); j++) {
-//			out << ", ";
-//			out << format(fmt,a[j]);
-//		}
-//		out << " ]";
-//		return out.str();
-//	}
-//}
-//
-//template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A) {
-//	return vec_to_string(A, "%7.3g");
-//}
-//
-//template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A, const char *fmt) {
-//	std::stringstream out;
-//	for (size_t j = 0; j < A.size(); j++) {
-//		out << vec_to_string(A[j], fmt);
-//    }
-//    return out.str();
-//}
-

 }; /* namespace CoolProp */
+
+
+
+
+#ifdef ENABLE_CATCH
+#include <math.h>
+#include <iostream>
+#include "catch.hpp"
+
+TEST_CASE("Internal consistency checks and example use cases for MatrixMath.h","[MatrixMath]")
+{
+	bool PRINT = false;
+
+	/// Test case for "SylthermXLT" by "Dow Chemicals"
+	std::vector<double> cHeat;
+	cHeat.clear();
+	cHeat.push_back(+1.1562261074E+03);
+	cHeat.push_back(+2.0994549103E+00);
+	cHeat.push_back(+7.7175381057E-07);
+	cHeat.push_back(-3.7008444051E-20);
+
+	std::vector<std::vector<double> > cHeat2D;
+	cHeat2D.push_back(cHeat);
+	cHeat2D.push_back(cHeat);
+
+	SECTION("Pretty printing tests") {
+
+		Eigen::MatrixXd matrix = Eigen::MatrixXd::Random(4,1);
+		std::string tmpStr;
+		if (PRINT) std::cout << std::endl;
+
+		CHECK_NOTHROW( tmpStr = CoolProp::vec_to_string(cHeat[0]) );
+		if (PRINT) std::cout << tmpStr << std::endl;
+		CHECK_NOTHROW( tmpStr = CoolProp::vec_to_string(cHeat) );
+		if (PRINT) std::cout << tmpStr << std::endl;
+		CHECK_NOTHROW( tmpStr = CoolProp::vec_to_string(cHeat2D) );
+		if (PRINT) std::cout << tmpStr << std::endl;
+
+		CHECK_NOTHROW( tmpStr = CoolProp::mat_to_string(CoolProp::vec_to_eigen(cHeat[0])) );
+		if (PRINT) std::cout << tmpStr << std::endl;
+		CHECK_NOTHROW( tmpStr = CoolProp::mat_to_string(CoolProp::vec_to_eigen(cHeat, 1)) );
+		if (PRINT) std::cout << tmpStr << std::endl;
+		CHECK_NOTHROW( tmpStr = CoolProp::mat_to_string(CoolProp::vec_to_eigen(cHeat, 2)) );
+		if (PRINT) std::cout << tmpStr << std::endl;
+		CHECK_NOTHROW( tmpStr = CoolProp::mat_to_string(CoolProp::vec_to_eigen(cHeat2D)) );
+		if (PRINT) std::cout << tmpStr << std::endl;
+	}
+
+	SECTION("Matrix modifications") {
+		Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(cHeat2D);
+
+		std::string tmpStr;
+		std::vector<std::vector<double> > vec2D;
+		if (PRINT) std::cout << std::endl;
+
+		CHECK_NOTHROW( CoolProp::removeColumn(matrix,1) );
+		if (PRINT) std::cout << CoolProp::mat_to_string(matrix) << std::endl;
+
+		CHECK_NOTHROW( CoolProp::removeRow(matrix,1) );
+		if (PRINT) std::cout << CoolProp::mat_to_string(matrix) << std::endl;
+
+		CHECK_THROWS( CoolProp::removeColumn(matrix,10) );
+		CHECK_THROWS( CoolProp::removeRow(matrix,10) );
+	}
+
+	SECTION("std::vector to Eigen::Matrix and back") {
+		std::vector<std::vector<double> > vec2D(cHeat2D);
+		Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(vec2D);
+		for (size_t i = 0; i < matrix.cols(); ++i) {
+			for (size_t j = 0; j < matrix.rows(); ++j) {
+				CHECK( fabs(matrix(j,i)-vec2D[j][i]) <= 1e-10 );
+			}
+		}
+		vec2D = CoolProp::eigen_to_vec(matrix);
+		for (size_t i = 0; i < matrix.cols(); ++i) {
+			for (size_t j = 0; j < matrix.rows(); ++j) {
+				CHECK( fabs(matrix(j,i)-vec2D[j][i]) <= 1e-10 );
+			}
+		}
+		std::vector<double> vec1D(cHeat);
+		matrix = CoolProp::vec_to_eigen(vec1D);
+		for (size_t i = 0; i < matrix.cols(); ++i) {
+			for (size_t j = 0; j < matrix.rows(); ++j) {
+				CHECK( fabs(matrix(j,i)-vec1D[j]) <= 1e-10 );
+			}
+		}
+		vec1D = CoolProp::eigen_to_vec1D(matrix);
+		for (size_t i = 0; i < matrix.cols(); ++i) {
+			for (size_t j = 0; j < matrix.rows(); ++j) {
+				CHECK( fabs(matrix(j,i)-vec1D[j]) <= 1e-10 );
+			}
+		}
+	}
+}
+
+
+#endif /* ENABLE_CATCH */
+
+
+//#include <unsupported/Eigen/Polynomials>
+//#include <iostream>
+////using namespace Eigen;
+////using namespace std;
+//
+//#include <vector>
+//#include <string>
+//#include <MatrixMath.h>
+//
+//int main()
+//{
+//Eigen::Vector4d roots = Eigen::Vector4d::Random();
+//std::cout << "Roots: " << roots.transpose() << std::endl;
+//Eigen::Matrix<double,5,1> polynomial;
+//Eigen::roots_to_monicPolynomial( roots, polynomial );
+//std::cout << "Polynomial: ";
+//for( int i=0; i<4; ++i ){ std::cout << polynomial[i] << ".x^" << i << "+ "; }
+//std::cout << polynomial[4] << ".x^4" << std::endl;
+//Eigen::Vector4d evaluation;
+//for( int i=0; i<4; ++i ){
+//evaluation[i] = Eigen::poly_eval( polynomial, roots[i] ); }
+//std::cout << "Evaluation of the polynomial at the roots: " << evaluation.transpose() << std::endl;
+//std::cout << std::endl;
+////
+////Eigen::MatrixXd coeffs = Eigen::MatrixXd::Random(5,1);
+////Eigen::MatrixXd input  = Eigen::MatrixXd::Random(2,1)*1e0;
+//Eigen::Vector4d coeffs = Eigen::Vector4d::Random()*1e2;
+//double input  = 1.9e0;
+//std::cout << "Coeffs: " << std::endl << coeffs.transpose() << std::endl;
+//double eval = Eigen::poly_eval( coeffs, input);
+//std::cout << "Evaluation of the polynomial at " << input << std::endl;
+//std::cout << eval << std::endl;
+//
+//double vec0 = 0.1;
+//std::vector<double> vec1(2,0.44);
+//std::vector< std::vector<double> > vec2;
+//vec2.push_back(std::vector<double>(2,0.2));
+//vec2.push_back(std::vector<double>(2,0.3));
+//
+//std::cout << CoolProp::vec_to_string(vec0) << std::endl;
+//std::cout << CoolProp::vec_to_string(vec1) << std::endl;
+//std::cout << CoolProp::vec_to_string(vec2) << std::endl;
+//
+//Eigen::Matrix<double,2,2> mat;
+//mat.setConstant(2,2,0.25);
+//std::vector< std::vector<double> > vec;
+//
+//CoolProp::convert(mat, vec);
+//std::cout << CoolProp::vec_to_string(vec) << std::endl;
+//
+////Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> mat;
+////mat.resize(6,2);
+//
+//Eigen::Matrix<double,2,2> mat2;
+//CoolProp::convert(vec2, mat2);
+//CoolProp::convert(mat2, vec);
+//std::cout << CoolProp::vec_to_string(vec) << std::endl;
+//
+//Eigen::Matrix<double,2,1> mat1;
+//CoolProp::convert(vec1, mat1);
+//std::vector<double> vec3;
+//CoolProp::convert(mat1, vec);
+//std::cout << CoolProp::vec_to_string(vec) << std::endl;
+//
+////std::vector< std::vector<double> > vec(vec2);
+////CoolProp::convert(mat,vec);
+//
+////std::cout << CoolProp::vec_to_string() << std::endl;
+//
+////Eigen::Matrix2d mat2 = CoolProp::convert(vec2);
+//
+////Eigen::MatrixXd mat2(10,10);
+////CoolProp::convert(vec2, mat2);
+//
+////std::cout << CoolProp::vec_to_string(CoolProp::convert(mat2)) << std::endl;
+//
+//}
--- a/src/Solvers.cpp
+++ b/src/Solvers.cpp
@@ -72,8 +72,13 @@ double Newton(FuncWrapper1D &f, double x0, double ftol, int maxiter, std::string
    {
 		fval = f.call(x);
 		dx = -fval/f.deriv(x);
-		x += dx;
-		
+
+		if (!ValidNumber(fval)){
+            throw ValueError("Residual function in newton returned invalid number");
+        };
+
+        x += dx;
+
 		if (fabs(dx/x) < 10*DBL_EPSILON)
 		{
 			return x;
@@ -102,6 +107,11 @@ In the secant function, a 1-D Newton-Raphson solver is implemented.  An initial
 */
 double Secant(FuncWrapper1D &f, double x0, double dx, double tol, int maxiter, std::string &errstring)
 {
+    #if defined(COOLPROP_DEEP_DEBUG)
+    static std::vector<double> xlog, flog;
+    xlog.clear(); flog.clear();
+    #endif
+
 	double x1=0,x2=0,x3=0,y1=0,y2=0,x,fval=999;
    int iter=1;
 	errstring = "";
@@ -112,8 +122,17 @@ double Secant(FuncWrapper1D &f, double x0, double dx, double tol, int maxiter, s
        if (iter==1){x1=x0; x=x1;}
        if (iter==2){x2=x0+dx; x=x2;}
        if (iter>2) {x=x2;}
-			fval=f.call(x);
-            if (!ValidNumber(fval)){throw ValueError("Residual function in secant returned invalid number");};
+
+            fval = f.call(x);
+
+            #if defined(COOLPROP_DEEP_DEBUG)
+                xlog.push_back(x);
+                flog.push_back(fval);
+            #endif
+
+            if (!ValidNumber(fval)){
+                throw ValueError("Residual function in secant returned invalid number");
+            };
        if (iter==1){y1=fval;}
        if (iter>1)
        {
--- a/src/Tests/CoolProp-Tests.cpp
+++ b/src/Tests/CoolProp-Tests.cpp
@@ -219,7 +219,7 @@ class TransportValidationFixture
 {
 protected:
    long double actual, x1, x2;
-    std::tr1::shared_ptr<CoolProp::AbstractState> pState;
+    shared_ptr<CoolProp::AbstractState> pState;
    int pair;
 public:
    TransportValidationFixture(){ }
@@ -504,7 +504,7 @@ class ConsistencyFixture
 {
 protected:
    long double hmolar, pmolar, smolar, umolar, rhomolar, T, p, x1, x2;
-    std::tr1::shared_ptr<CoolProp::AbstractState> pState;
+    shared_ptr<CoolProp::AbstractState> pState;
    int pair;
 public:
    ConsistencyFixture(){}
--- a/src/Tests/test_main.cxx
+++ b/src/Tests/test_main.cxx
@@ -0,0 +1,5 @@
+#ifdef ENABLE_CATCH
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
+#endif /* ENABLE_CATCH */
--- a/src/main.cxx
+++ b/src/main.cxx
@@ -18,6 +18,8 @@ using namespace CoolProp;
 #endif
 #include "SpeedTest.h"

+#include "crossplatform_shared_ptr.h"
+
 //#include <vld.h>

 void generate_melting_curve_data(const char* file_name, const char *fluid_name, double Tmin, double Tmax)
@@ -25,7 +27,7 @@ void generate_melting_curve_data(const char* file_name, const char *fluid_name,

    FILE *fp;
    fp = fopen(file_name,"w");
-    std::tr1::shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"),std::string(fluid_name)));
+    shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"),std::string(fluid_name)));
    for (double T = Tmin; T < Tmax; T += 0.1)
    {
        try{
@@ -55,9 +57,74 @@ struct element
            double d,t,ld;
            int l;
        };
+
 int main()
 {
-    set_debug_level(0);
+    set_debug_level(1);
+
+
+    if (0)
+    {
+        // First type (slowest, most string processing, exposed in DLL)
+        double r0A = PropsSI("Dmolar","T",298,"P",1e5,"Propane[0.5]&Ethane[0.5]"); // Default backend is HEOS
+        double r0B = PropsSI("Dmolar","T",298,"P",1e5,"HEOS::Propane[0.5]&Ethane[0.5]");
+        double r0C = PropsSI("Dmolar","T",298,"P",1e5,"REFPROP::Propane[0.5]&Ethane[0.5]");
+
+        std::vector<double> z(2,0.5);
+        // Second type (C++ only, a bit faster)
+        double r1A = PropsSI("Dmolar","T",298,"P",1e5,"Propane&Ethane", z);
+        double r1B = PropsSI("Dmolar","T",298,"P",1e5,"HEOS::Propane&Ethane", z);
+        double r1C = PropsSI("Dmolar","T",298,"P",1e5,"REFPROP::Propane&Ethane", z);
+
+        //const double *pz = &(z[0]);
+        //int n = z.size();
+        //// Third type (DLL)
+        //double r2A = PropsSIZ("Dmolar","T",298,"P",1e5,"Propane&Ethane", pz, n);
+        //double r2B = PropsSIZ("Dmolar","T",298,"P",1e5,"HEOS::Propane&Ethane", pz, n);
+        //double r2C = PropsSIZ("Dmolar","T",298,"P",1e5,"REFPROP::Propane&Ethane", pz, n);
+
+        double tt = 0;
+    }
+
+    if (0)
+    {
+        shared_ptr<CoolProp::AbstractState> AS(AbstractState::factory("HEOS","ETHANE&PROPANE"));
+        std::vector<double>x(2,0.5);
+        AS->set_mole_fractions(x);
+        //AS->build_phase_envelope("");
+        AS->update(PQ_INPUTS,648000,0);
+        for (int Q = 0; Q <= 1; Q++)
+        {
+            std::vector<double> TT,PP,RR;
+            for (double p = 101325; ;p*=1.00005)
+            {
+                try{
+                AS->update(PQ_INPUTS,p,Q);
+                double T = AS->T();
+                double rho = AS->rhomolar();
+                double pp = AS->p();
+                TT.push_back(T);
+                PP.push_back(pp);
+                RR.push_back(rho);
+                printf("%g, %g, %g\n", T, rho, pp);
+                }
+                catch(std::exception &e){std::cout << e.what() << std::endl; break;}
+                catch(...){break;}
+            }
+            rapidjson::Document d;
+            d.SetObject();
+            cpjson::set_double_array("T",TT,d,d);
+            cpjson::set_double_array("P",PP,d,d);
+            cpjson::set_double_array("rho",RR,d,d);
+            std::string fname = format("%d.json",Q);
+            FILE *fp = fopen(fname.c_str(),"w");
+            fprintf(fp,"%s",cpjson::json2string(d).c_str());
+            fclose(fp);
+        }
+
+        double rr = 0;
+        return 0;
+    }
    if (0)
    {
        std::string NBP_refs[] = {"D5","D6","MD2M","MDM","Benzene","Helium","Ethylene","Ethanol","n-Dodecane","Benzene","n-Undecane","Neon","Fluorine","Methanol","Acetone","Methane","Ethane","n-Pentane","n-Hexane","n-Heptane","n-Octane","CycloHexane","MD3M","MM","D4","MethylPalmitate","MethylStearate","MethylOleate","MethylLinoleate","MethylLinolenate","m-Xylene","Air"};
@@ -154,24 +221,46 @@ int main()
            std::cout << format("%s %17.15g\n", S->name().c_str(), S->p());
        }
    }
+    if (1){
+        shared_ptr<CoolProp::AbstractState> ASR(CoolProp::AbstractState::factory("HEOS","H2S"));
+        ASR->update(PT_INPUTS, 1000e6, 200);
+        double v  = ASR->viscosity();
+        int rr =0;
+    }
    if (0)
    {
-        double rrr0 = PropsSI("P","T",200,"Dmolar",14000,"REFPROP::R125");
-        double rrr2 = PropsSI("P","T",200,"Dmolar",14000,"R125");
+        shared_ptr<CoolProp::AbstractState> ASR(CoolProp::AbstractState::factory("REFPROP","CO2"));
+        double p1 = ASR->calc_melt_p_T(250);
+
+        shared_ptr<CoolProp::AbstractState> ASC(CoolProp::AbstractState::factory("HEOS","CO2"));
+        double p2 = ASC->calc_melt_p_T(250);
+
+        double rrr1 = PropsSI("D","T",200,"P",300,"Water");
+        double rrr0 = PropsSI("Cvmolar","T",200,"Dmolar",14000,"REFPROP::R125");
+        double rrr2 = PropsSI("speed_of_sound","T",300,"Dmolar",700,"R125");
        double rrr =0 ;
    }
    if (0)
    {
-        std::tr1::shared_ptr<AbstractState> ASR(AbstractState::factory("REFPROP","CO2"));
+        shared_ptr<AbstractState> ASR(AbstractState::factory("REFPROP","CO2"));
        ASR->update(QT_INPUTS, 1, 304);
        double muR0 = ASR->conductivity();

-        std::tr1::shared_ptr<AbstractState> ASC(AbstractState::factory("HEOS","CO2"));
+        shared_ptr<AbstractState> ASC(AbstractState::factory("HEOS","CO2"));
        ASC->update(QT_INPUTS, 1, 304);
        double muC = ASC->conductivity();
        double rr = 4;
    }
-    if (1)
+    if (0)
+    {
+        #if ENABLE_CATCH
+            std::vector<std::string> tags;
+            tags.push_back("[transport]");
+            run_user_defined_tests(tags);
+            double rr = 0;
+        #endif
+    }
+    if (0)
    {
        /*double h1 = PropsSI("S","P",101325,"Q",0,"n-Pentane");
        std::string er = get_global_param_string("errstring");
@@ -181,18 +270,15 @@ int main()
        //std::string RPname = get_fluid_param_string("Water", "REFPROPname");
        //std::string s = get_BibTeXKey("n-Propane", "rr");

-        double rr0 = PropsSI("L","T",647.35,"Dmass",272,"Water");
+

        #if ENABLE_CATCH
-        std::vector<std::string> tags;
-        tags.push_back("[rp1485]");
-        //run_user_defined_tests(tags);
        run_tests();
        #endif

        std::string fl = get_global_param_string("FluidsList");
        double rr = PropsSI("D", "P", 3e5, "T", 300, "Nitrogen");
-        std::tr1::shared_ptr<AbstractState> AS(AbstractState::factory("HEOS","Nitrogen"));
+        shared_ptr<AbstractState> AS(AbstractState::factory("HEOS","Nitrogen"));

        AS->update(DmolarT_INPUTS, 40, 300);
        double p1 = AS->umolar();
@@ -246,11 +332,13 @@ int main()
    if (0)
    {

-        double n[] = {0.0125335479355233,                        7.8957634722828,                        -8.7803203303561,                        0.31802509345418,                        -0.26145533859358,                        -0.0078199751687981,                        0.0088089493102134,                        -0.66856572307965,                        0.20433810950965,                        -6.621260503968699e-005,                        -0.19232721156002,                        -0.25709043003438,                        0.16074868486251,                        -0.040092828925807,                        3.9343422603254e-007,                        -7.5941377088144e-006,                        0.00056250979351888,                        -1.5608652257135e-005,                        1.1537996422951e-009,                        3.6582165144204e-007,                        -1.3251180074668e-012,                        -6.2639586912454e-010,                        -0.10793600908932,                        0.017611491008752,                        0.22132295167546,                        -0.40247669763528,                        0.58083399985759,                        0.0049969146990806,                        -0.031358700712549,                        -0.74315929710341,                        0.4780732991548,                        0.020527940895948,                        -0.13636435110343,                        0.014180634400617,                        0.008332650488071301,                        -0.029052336009585,                        0.038615085574206,                        -0.020393486513704,                        -0.0016554050063734,                        0.0019955571979541,                        0.00015870308324157,                        -1.638856834253e-005,                        0.043613615723811,                        0.034994005463765,                        -0.076788197844621,                        0.022446277332006,                        -6.2689710414685e-005,                        -5.5711118565645e-010,                        -0.19905718354408,                        0.31777497330738,                        -0.11841182425981  };
-        double d[] = {1,                         1,                        1,                        2,                        2,                        3,                        4,                        1,                        1,                        1,                        2,                        2,                        3,                        4,                        4,                        5,                        7,                        9,                        10,                        11,                        13,                        15,                        1,                        2,                        2,                        2,                        3,                        4,                        4,                        4,                        5,                        6,                        6,                        7,                        9,                        9,                        9,                        9,                        9,                        10,                        10,                        12,                        3,                        4,                        4,                        5,                        14,                        3,                        6,                        6,                        6                    };
-        double t[] = {-0.5,                        0.875,                        1,                        0.5,                        0.75,                        0.375,                        1,                        4,                        6,                        12,                        1,                        5,                        4,                        2,                        13,                        9,                        3,                        4,                        11,                        4,                        13,                        1,                        7,                        1,                        9,                        10,                        10,                        3,                        7,                        10,                        10,                        6,                        10,                        10,                        1,                        2,                        3,                        4,                        8,                        6,                        9,                        8,                        16,                        22,                        23,                        23,                        10,                        50,                        44,                        46,                        50                    };
-        double l[] = {0,                        0,                        0,                        0,                        0,                        0,                        0,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        3,                        3,                        3,                        3,                        4,                        6,                        6,                        6,                        6                    };
-        double summer = 0;
+
+        typedef double dbltype;
+        dbltype n[] = {0.0125335479355233,                        7.8957634722828,                        -8.7803203303561,                        0.31802509345418,                        -0.26145533859358,-0.0078199751687981,0.0088089493102134,                        -0.66856572307965,                        0.20433810950965,                        -6.621260503968699e-005,                        -0.19232721156002,                        -0.25709043003438,                        0.16074868486251,                        -0.040092828925807,                        3.9343422603254e-007,                        -7.5941377088144e-006,                        0.00056250979351888,                        -1.5608652257135e-005,                        1.1537996422951e-009,                        3.6582165144204e-007,                        -1.3251180074668e-012,                        -6.2639586912454e-010,                        -0.10793600908932,                        0.017611491008752,                        0.22132295167546,                        -0.40247669763528,                        0.58083399985759,                        0.0049969146990806,                        -0.031358700712549,                        -0.74315929710341,                        0.4780732991548,                        0.020527940895948,                        -0.13636435110343,                        0.014180634400617,                        0.008332650488071301,                        -0.029052336009585,                        0.038615085574206,                        -0.020393486513704,                        -0.0016554050063734,                        0.0019955571979541,                        0.00015870308324157,                        -1.638856834253e-005,                        0.043613615723811,                        0.034994005463765,-0.076788197844621,0.022446277332006,-6.2689710414685e-005,-5.5711118565645e-010,-0.19905718354408,0.31777497330738,-0.11841182425981};
+        dbltype d[] = {1,                         1,                        1,                        2,                        2,                        3,                        4,                        1,                        1,                        1,                        2,                        2,                        3,                        4,                        4,                        5,                        7,                        9,                        10,                        11,                        13,                        15,                        1,                        2,                        2,                        2,                        3,                        4,                        4,                        4,                        5,                        6,                        6,                        7,                        9,                        9,                        9,                        9,                        9,                        10,                        10,                        12,                        3,                        4,                        4,                        5,                        14,                        3,                        6,                        6,                        6                    };
+        dbltype t[] = {-0.5,                        0.875,                        1,                        0.5,                        0.75,                        0.375,                        1,                        4,                        6,                        12,                        1,                        5,                        4,                        2,                        13,                        9,                        3,                        4,                        11,                        4,                        13,                        1,                        7,                        1,                        9,                        10,                        10,                        3,                        7,                        10,                        10,                        6,                        10,                        10,                        1,                        2,                        3,                        4,                        8,                        6,                        9,                        8,                        16,                        22,                        23,                        23,                        10,                        50,                        44,                        46,                        50                    };
+        dbltype l[] = {0,                        0,                        0,                        0,                        0,                        0,                        0,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        1,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        2,                        3,                        3,                        3,                        3,                        4,                        6,                        6,                        6,                        6                    };
+        dbltype summer = 0;
        std::vector<element> elements;
        for (std::size_t i = 0; i < 51; ++i)
        {
@@ -258,42 +346,54 @@ int main()
            el.d = d[i];
            el.t = t[i];
            el.l = (int)l[i];
-            el.ld = (double)l[i];
+            el.ld = (dbltype)l[i];
            elements.push_back(el);
        }

        long N = 1000000;
+        dbltype pow_delta_li, di, ti, lid;
+        std::vector<dbltype> s(51);
+        dbltype delta, log_tau, log_delta, tau, gamma;
        double t1 = clock();
+
        for (std::size_t ii = 0; ii < N; ++ii)
        {
-            double delta = 1.3, tau = 0.7;
-            double log_tau  = log(tau), log_delta = log(delta);
+            delta = 1.3, tau = 0.7;
+            log_tau  = log(tau), log_delta = log(delta);
+            pow_delta_li = pow(delta, 5);
            for (int jj = 0; jj < 51; ++jj)
            {
-                double di = d[jj], ti = t[jj], lid = l[jj];
-                int li = (int)lid;
-                double pow_delta_li;
-                if (li > 0){
-                    pow_delta_li = pow(delta, li);
-                    summer += (di-lid*pow_delta_li)*exp(ti*log_tau+(di-1)*log_delta-pow_delta_li);
-                }
-                else{
-                    summer += di*exp(ti*log_tau+(di-1)*log_delta);
-                }
+                //di = d[jj]; ti = t[jj]; lid = l[jj]; gamma = 1;
+                //int li = (int)lid;
+
+                element &el = elements[jj];
+                ti = el.t;
+                //int li = (int)lid;
+
+                summer += exp(ti*log_tau);//(di-gamma*lid*pow_delta_li)*pow(tau,ti)*pow_delta_li*exp(-gamma*pow_delta_li);
+                //summer += __ieee754_exp(ti*log_tau);//(di-gamma*lid*pow_delta_li)*pow(tau,ti)*pow_delta_li*exp(-gamma*pow_delta_li);
+//                if (li > 0){
+//
+//
+//                }
+//                else{
+//                    s[jj] += di*exp(ti*log_tau+(di-1)*log_delta);
+//                }
            }
        }
+        //summer = std::accumulate(s.begin(), s.end(), (dbltype)(0));
        double t2 = clock();
-        double elap = (t2-t1)/CLOCKS_PER_SEC/((double)N)*1e6;
+        double elap = (t2-t1)/CLOCKS_PER_SEC/((dbltype)N)*1e6;
        printf("%g %g\n",elap, summer);
-
+        int rr =5;
    }
-    if (0)
+    if (1)
    {
-        std::tr1::shared_ptr<AbstractState> MixRP(AbstractState::factory(std::string("REFPROP"),std::string("propane")));
+        shared_ptr<AbstractState> MixRP(AbstractState::factory(std::string("REFPROP"),std::string("propane")));
        MixRP->update(QT_INPUTS, 0, 330);
        long double s1 = MixRP->surface_tension();

-        std::tr1::shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("propane")));
+        shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("propane")));
        Mix->update(QT_INPUTS, 0, 330);
        long double s2 = Mix->surface_tension();
    }
@@ -302,7 +402,7 @@ int main()

        double T = 300;

-        std::tr1::shared_ptr<AbstractState> MixRP(AbstractState::factory(std::string("REFPROP"),std::string("propane")));
+        shared_ptr<AbstractState> MixRP(AbstractState::factory(std::string("REFPROP"),std::string("propane")));
        {
            long N = 100000;
            double t1 = clock(), summer = 0;
@@ -320,7 +420,7 @@ int main()
        double cp2 = MixRP->cpmolar();
        double T2 = MixRP->T();

-        std::tr1::shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("propane")));
+        shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("propane")));
        {
            long N = 100000;
            double t1 = clock(), summer = 0;
@@ -348,7 +448,7 @@ int main()

        int inputs = PQ_INPUTS; double val1 = p, val2 = Q;

-        std::tr1::shared_ptr<AbstractState> MixRP(AbstractState::factory(std::string("REFPROP"), std::string("Ethane,propane")));
+        shared_ptr<AbstractState> MixRP(AbstractState::factory(std::string("REFPROP"), std::string("Ethane,propane")));
        MixRP->set_mole_fractions(z);
        MixRP->update(inputs, val1, val2);
        double p2 = MixRP->p();
@@ -361,7 +461,7 @@ int main()
        double phi20 = MixRP->fugacity_coefficient(0);
        double phi21 = MixRP->fugacity_coefficient(1);

-        std::tr1::shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("Ethane,propane")));
+        shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("Ethane,propane")));
        Mix->set_mole_fractions(z);
        Mix->update(inputs, val1, val2);
        double p1 = Mix->p();
@@ -381,7 +481,7 @@ int main()
        std::vector<long double> z(N, 1.0/N);
        double Q = 0, T = 250, p = 300000;

-        std::tr1::shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("Ethane,n-Propane")));
+        shared_ptr<AbstractState> Mix(AbstractState::factory(std::string("HEOS"), std::string("Ethane,n-Propane")));

        Mix->set_mole_fractions(z);

@@ -396,7 +496,7 @@ int main()
        time_t t1,t2;

        std::size_t N = 1000000;
-        std::tr1::shared_ptr<AbstractState> State(AbstractState::factory(std::string("HEOS"), std::string("Water")));
+        shared_ptr<AbstractState> State(AbstractState::factory(std::string("HEOS"), std::string("Water")));

        double p = State->p();
        double summer = 0;
@@ -425,7 +525,7 @@ int main()

    if (0)
    {
-        std::tr1::shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"), std::string("Methane|Ethane")));
+        shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"), std::string("Methane|Ethane")));

        std::vector<long double> x(2,0.5);
        State->set_mole_fractions(x);
@@ -441,7 +541,7 @@ int main()
        long N = 100000;
        for (long ii = 0; ii < N; ii++)
        {
-            std::tr1::shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"), std::string("Methane")));
+            shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"), std::string("Methane")));
            //AbstractState *State = new REFPROPBackend("Methane");
        }
        t2 = clock();
@@ -451,7 +551,7 @@ int main()

    if(0)
    {
-        std::tr1::shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"), std::string("Methane")));
+        shared_ptr<AbstractState> State(AbstractState::factory(std::string("REFPROP"), std::string("Methane")));

        State->update(DmassT_INPUTS,1,300);
        double hh = State->hmolar();
--- a/wrappers/C#/Example.cs
+++ b/wrappers/C#/Example.cs
@@ -1,6 +1,5 @@
 using System;
 using System.Collections.Generic;
-using System.Linq;
 using System.Text; 

 namespace ConsoleApplication1
--- a/wrappers/Octave/Example.m
+++ b/wrappers/Octave/Example.m
--- a/wrappers/Python/CoolProp5/AbstractState.pxd
+++ b/wrappers/Python/CoolProp5/AbstractState.pxd
@@ -20,11 +20,11 @@ cdef class AbstractState:
    cpdef double cpmolar(self)
    cpdef double cvmolar(self)
    cpdef double rhomass(self)
-    cpdef double hmass(self)
-    cpdef double smass(self)
-    cpdef double cpmass(self)
-    cpdef double cvmass(self)
-    cpdef double speed_sound(self)
+    cpdef double hmass(self) except *
+    cpdef double smass(self) except *
+    cpdef double cpmass(self) except *
+    cpdef double cvmass(self) except *
+    cpdef double speed_sound(self) except *
    
-    
-    cpdef double keyed_output(self, long)
+    cpdef double molar_mass(self) except *
+    cpdef double keyed_output(self, long) except *
--- a/wrappers/Python/CoolProp5/AbstractState.pyx
+++ b/wrappers/Python/CoolProp5/AbstractState.pyx
@@ -64,4 +64,6 @@ cdef class AbstractState:
        """ Get the speed of sound in m/s - wrapper of c++ function :cpapi:`AbstractState::speed_sound` """
        return self.thisptr.speed_sound()
    
-        
+    cpdef double molar_mass(self) except *: 
+        """ Get the molar mass kg/mol - wrapper of c++ function :cpapi:`AbstractState::molar_mass` """
+        return self.thisptr.molar_mass()    
--- a/wrappers/Python/CoolProp5/cAbstractState.pxd
+++ b/wrappers/Python/CoolProp5/cAbstractState.pxd
@@ -32,7 +32,13 @@ cdef extern from "AbstractState.h" namespace "CoolProp":
        double cvmass() except +
        
        double keyed_output(long) except+
+        double molar_mass() except+
+        double gas_constant() except+
+        double build_phase_envelope() except+
+        double viscosity() except+
+        double conductivity() except+
+        double surface_tension() except+

 # The static factory method for the AbstractState
 cdef extern from "AbstractState.h" namespace "CoolProp::AbstractState":
-    AbstractState* factory(const string &backend, const string &fluid_string)
+    AbstractState* factory(const string &backend, const string &fluid_string) except+
--- a/wrappers/StaticLibrary/readme.txt
+++ b/wrappers/StaticLibrary/readme.txt
@@ -1 +0,0 @@
-Sorry, all build scripts have been moved to the SharedLibrary folder. Please have a look there.
				`@@ -1 +0,0 @@`
				`Sorry, all build scripts have been moved to the SharedLibrary folder. Please have a look there.`