Added more fitting stuff, still some problems with the enumerations in C++

dev/incompressible_liquids/.gitignore

@@ -0,0 +1,2 @@
+/report/
+/.ipynb_checkpoints/

dev/incompressible_liquids/CPIncomp/BaseObjects.py

@@ -26,7 +26,10 @@ class IncompressibleData(object):
     SOURCE_COEFFS   = 'coefficients'
     SOURCE_NOT_SET  = 'notdefined'
     
-    maxLog = np.log(np.finfo(np.float64).max-1)
+    maxLin = np.finfo(np.float64).max-1
+    minLin = -maxLin
+    
+    maxLog = np.log(maxLin)
     minLog = -maxLog
     
     def __init__(self):                
@@ -36,8 +39,8 @@ class IncompressibleData(object):
         self.data   = None # None #np.zeros((10,10))
         self.xData  = None # In case you need a customised first data set (temperature?)
         self.yData  = None # In case you need a customised second data set (concentration?)
-        self.rsq    = None # Coefficient of determination
-        self.DEBUG = False
+        self.sErr   = None # Coefficient of determination
+        self.DEBUG  = False
         
     @staticmethod
     def baseFunc(x, y=0.0, xbase=0.0, ybase=0.0, eqnType=None, c=None):
@@ -77,10 +80,10 @@ class IncompressibleData(object):
     @staticmethod
     def baseLogexponential(co, x):
         r,c,coeffs = IncompressibleFitter.shapeArray(co)
-        if not ( (r==4 and c==1) or (r==1 and c==4) ):
+        if not ( (r==3 and c==1) or (r==1 and c==3) ):
             raise ValueError("You have to provide a 4,1 matrix of coefficients, not ({0},{1}).".format(r,c))
         coeffs_tmp = np.array(coeffs.flat)
-        return np.exp(np.clip(np.log(1.0/(x+coeffs_tmp[0]) + 1.0/(x+coeffs_tmp[1])**2)*coeffs_tmp[2]+coeffs_tmp[3],IncompressibleData.minLog,IncompressibleData.maxLog))
+        return np.exp(np.clip(np.log(np.clip(1.0/(x+coeffs_tmp[0]) + 1.0/(x+coeffs_tmp[0])**2,1e-10,IncompressibleData.maxLin))*coeffs_tmp[1]+coeffs_tmp[2],IncompressibleData.minLog,IncompressibleData.maxLog))
 
     @staticmethod
     def basePolyOffset(co, x):
@@ -110,40 +113,69 @@ class IncompressibleData(object):
         #res = None
         #r2 = None
         
-        res,r2 = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
+        res,sErr = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
                   xbase=xbase, ybase=ybase, \
                   eqnType=self.type, \
                   coeffs=self.coeffs, DEBUG=self.DEBUG)
         
-        #count = 0
-        #while r2<0.9 and count<1:
-            ##if self.DEBUG: print("Poor solution found, trying once more with more coefficients.")
-            #if self.type==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL:
-                #if self.DEBUG: print("Poor solution found, trying once more with log exponential.")
-                #self.type=IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL
-                #self.coeffs = np.array([-250,-250,2,2])
-                #res,r2 = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
-                      #xbase=xbase, ybase=ybase, \
-                      #eqnType=self.type, \
-                      #coeffs=self.coeffs, DEBUG=self.DEBUG)
+        bestCoeffs = np.copy(res)
+        bestType   = self.type
+        bestsErr   = np.copy(sErr)
+        bestRMS    = np.sqrt(np.square(bestsErr).mean()).sum()
+        
+        count = 0
+        while bestRMS>0.01 and count<2:
+            #if self.DEBUG: print("Poor solution found, trying once more with more coefficients.")
+            if self.type==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL:
+                if self.DEBUG: print("Poor solution found with exponential, trying once more with log exponential.")
+                self.type=IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL
+                self.coeffs = np.array([-250,1.5,10])
+                res,sErr = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
+                      xbase=xbase, ybase=ybase, \
+                      eqnType=self.type, \
+                      coeffs=self.coeffs, DEBUG=self.DEBUG)
             
-            #if self.type==IncompressibleData.INCOMPRESSIBLE_POLYNOMIAL or \
-              #self.type==IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL:
-                #if self.DEBUG: print("Poor solution found, trying once more with more coefficients.")
-                #self.coeffs = np.zeros(np.round(np.array(self.coeffs.shape) * 1.5))
-                #res,r2 = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
-                      #xbase=xbase, ybase=ybase, \
-                      #eqnType=self.type, \
-                      #coeffs=self.coeffs, DEBUG=self.DEBUG)
-            #count += 1
+            elif self.type==IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL:
+                if self.DEBUG: print("Poor solution found with log exponential, trying once more with exponential polynomial.")
+                self.type=IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL
+                self.coeffs = np.zeros((4,6))
+                res,sErr = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
+                      xbase=xbase, ybase=ybase, \
+                      eqnType=self.type, \
+                      coeffs=self.coeffs, DEBUG=self.DEBUG)
+                
+#            elif self.type==IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL:
+#                if self.DEBUG: print("Poor solution found with exponential polynomial, trying once more with normal polynomial.")
+#                self.type=IncompressibleData.INCOMPRESSIBLE_POLYNOMIAL
+#                self.coeffs = np.zeros((4,6))
+#                res,sErr = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
+#                      xbase=xbase, ybase=ybase, \
+#                      eqnType=self.type, \
+#                      coeffs=self.coeffs, DEBUG=self.DEBUG)
+
+            RMS = np.sqrt(np.square(sErr).mean()).sum()
+            if RMS<bestRMS: # Better fit
+                bestCoeffs = np.copy(res)
+                bestType   = self.type
+                bestsErr   = np.copy(sErr)
+                bestRMS    = RMS
+            
+            count += 1
          
-        if res==None:
+        if bestCoeffs==None:
             if self.DEBUG: print("There was a fitting error, no solution found.")
-        elif IncompressibleFitter.allClose(res, self.coeffs):
+        elif IncompressibleFitter.allClose(bestCoeffs, self.coeffs):
             if self.DEBUG: print("Coefficients did not change.")
         else:
-            self.coeffs = res
-            self.rsq    = r2 
+            if self.DEBUG: print("Best fit for: {0}".format(bestType))
+            self.coeffs = bestCoeffs
+            self.type   = bestType 
+            self.sErr   = bestsErr
+            
+            #if self.DEBUG: print("Fitting statistics:")
+            #SSE = np.square(self.sErr).sum() # Sum of squares due to error
+            #SST = ((zData-zData.mean())**2).sum()
+            #R2  = 1-(ssErr/ssTot )
             
 
     def setxData(self, xData):
@@ -257,7 +289,7 @@ class IncompressibleFitter(object):
         
         if zr==1 and zc==1: # 
             if DEBUG: print("Data no set, we cannot fit the coefficients")
-            return None
+            return None,None
                 
         if (xc!=1): raise ValueError("The first input has to be a 2D array with one column.")
         if (yr!=1): raise ValueError("The second input has to be a 2D array with one row.")
@@ -298,10 +330,12 @@ class IncompressibleFitter(object):
             z_input = np.copy(z)
             if eqnType==IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL:
                 z_input = np.log(z_input)
-            coeffs,r2 = IncompressibleFitter.getCoeffs2d(x_input, y_input, z_input, cr-1, cc-1, DEBUG=DEBUG)
+            coeffs,sErr = IncompressibleFitter.getCoeffs2d(x_input, y_input, z_input, cr-1, cc-1, DEBUG=DEBUG)
             if DEBUG: print("Coefficients after fitting: \n{0}".format(coeffs))
-            if DEBUG: print("Coefficient of determination: {0}".format(r2))
-            return coeffs,r2
+            if DEBUG: print("Standard deviation: {0}".format(np.nanstd(sErr)))
+            if DEBUG: print("Sum of squared errors: {0}".format(np.square(sErr).sum()))
+            if DEBUG: print("Root mean squared errors: {0}".format(np.sqrt(np.square(sErr).mean()).sum()))
+            return coeffs,sErr
         
         
         # Select if 1D or 2D fitting
@@ -309,14 +343,16 @@ class IncompressibleFitter(object):
             if DEBUG: print("1D function detected.")
             if yc==1: 
                 if DEBUG: print("Fitting {0} in x-direction.".format(eqnType))
-                coeffs,r2 = IncompressibleFitter.getCoeffsIterative1D(x, z, eqnType=eqnType, coeffs=coeffs, DEBUG=DEBUG)
+                coeffs,sErr = IncompressibleFitter.getCoeffsIterative1D(x, z, eqnType=eqnType, coeffs=coeffs, DEBUG=DEBUG)
             elif xr==1: 
                 if DEBUG: print("Fitting {0} in y-direction.".format(eqnType))
-                coeffs,r2 = IncompressibleFitter.getCoeffsIterative1D(y, z, eqnType=eqnType, coeffs=coeffs, DEBUG=DEBUG)
+                coeffs,sErr = IncompressibleFitter.getCoeffsIterative1D(y, z, eqnType=eqnType, coeffs=coeffs, DEBUG=DEBUG)
             else: raise ValueError("Unknown error in matrix shapes.")    
             if DEBUG: print("Coefficients after fitting: \n{0}".format(coeffs))
-            if DEBUG: print("Coefficient of determination: {0}".format(r2))
-            return coeffs, r2
+            if DEBUG: print("Standard deviation:    {0}".format(np.nanstd(sErr)))
+            if DEBUG: print("Sum of squared errors: {0}".format(np.square(sErr).sum()))
+            if DEBUG: print("Root mean squared errors: {0}".format(np.sqrt(np.square(sErr).mean()).sum()))
+            return coeffs, sErr
 
         elif yc>1: # 2D fitting
             raise ValueError("There are no other 2D fitting functions than polynomials, cannot use {0}.".format(eqnType))
@@ -345,8 +381,6 @@ class IncompressibleFitter(object):
         if y_order==None: raise ValueError("You did not provide data for y_order.")
         if DEBUG==None: raise ValueError("You did not provide data for DEBUG.")
         
-        r2 = None
-        
         x_order += 1
         y_order += 1
         #To avoid overfitting, we only use the upper left triangle of the coefficient matrix
@@ -387,6 +421,8 @@ class IncompressibleFitter(object):
         # Remove np.nan elements
         mask = np.isfinite(zz)
         A = A[mask]
+        xx = xx[mask]
+        yy = yy[mask]
         zz = zz[mask]
         
         if (len(A) < cols):
@@ -399,10 +435,10 @@ class IncompressibleFitter(object):
         if DEBUG: print(rank)
         if DEBUG: print(singulars)
         
-        if resids.size>0:
-            r2 = 1 - resids / (zz.size * zz.var())
-        else:
-            r2 = 0
+        #if resids.size>0:
+        #    r2 = 1 - resids / (zz.size * zz.var())
+        #else:
+        #    r2 = 0
         #print("\n r2 2d: ",r2.shape,r2,"\n")
         
         #Rearrange coefficients to a matrix shape
@@ -410,7 +446,8 @@ class IncompressibleFitter(object):
         for i, (xi,yi) in enumerate(xy_exp): # makes columns
             C[xi][yi] = coeffs[i]
         
-        return C, r2
+        sErr = zz - np.polynomial.polynomial.polyval2d(xx, yy, C)
+        return C, sErr
     
     @staticmethod
     def getCoeffsIterative1D(x_in, z_in, eqnType, coeffs, DEBUG=False):
@@ -421,7 +458,7 @@ class IncompressibleFitter(object):
         if coeffs==None: raise ValueError("You did not provide data for the coefficients.")
         if DEBUG==None: raise ValueError("You did not provide data for DEBUG.")
         
-        r2 = None
+        sErr = None 
         
         #fit = "Powell" # use Powell's algorithm
         #fit = "BFGS" # use Broyden-Fletcher-Goldfarb-Shanno
@@ -483,16 +520,17 @@ class IncompressibleFitter(object):
                     if np.any(popt!=coeffs):
                         success = True
                         if DEBUG: print("Fit succeeded with: {0}".format(algorithm))
-#                        print "Fit succeeded for "+fit[counter]+": "
-#                        print "data: {0}, func: {1}".format(yData[ 2],func(xData[ 2], popt))
-#                        print "data: {0}, func: {1}".format(yData[ 6],func(xData[ 6], popt))
-#                        print "data: {0}, func: {1}".format(yData[-1],func(xData[-1], popt))
+                        sErr = zData - func(xData, popt)
+                        #print "Fit succeeded for "+fit[counter]+": "
+                        #print "data: {0}, func: {1}".format(yData[ 2],func(xData[ 2], popt))
+                        #print "data: {0}, func: {1}".format(yData[ 6],func(xData[ 6], popt))
+                        #print "data: {0}, func: {1}".format(yData[-1],func(xData[-1], popt))
                         #if DEBUG: print("Estimated covariance of parameters: {0}".format(pcov))
-                        ssErr = np.sqrt(np.diag(pcov)).sum()
-                        ssTot = ((zData-zData.mean())**2).sum()
-                        r2 = 1-(ssErr/ssTot )
+                        #ssErr = np.sqrt(np.diag(pcov)).sum()
+                        #ssTot = ((zData-zData.mean())**2).sum()
+                        #r2 = 1-(ssErr/ssTot )
                         #print("\n r2 FMA: ",r2.shape,r2,"\n")
-                        return popt,r2
+                        return popt,sErr
                     else:
                         if DEBUG: print("Fit failed for {0}.".format(algorithm))
                         if DEBUG: sys.stdout.flush()
@@ -514,12 +552,13 @@ class IncompressibleFitter(object):
                     if res.success:
                         success = True
                         if DEBUG: print("Fit succeeded with: {0}".format(algorithm))
-                        if res.has_key('fvec'):
-                            ssErr = (res['fvec']**2).sum()
-                            ssTot = ((zData-zData.mean())**2).sum()
-                            r2 = 1-(ssErr/ssTot )
+                        sErr = zData - fun(np.array(res.x), xData)
+                        #if res.has_key('fvec'):
+                            #ssErr = (res['fvec']**2).sum()
+                            #ssTot = ((zData-zData.mean())**2).sum()
+                            #r2 = 1-(ssErr/ssTot )
                         #print("\n r2 : ",r2.shape,r2,algorithm,"\n")                    
-                        return res.x,r2
+                        return res.x,sErr
                     else:
                         if DEBUG: print("Fit failed for {0}.".format(algorithm))
                         if DEBUG: sys.stdout.flush()
@@ -546,7 +585,7 @@ class IncompressibleFitter(object):
                 if DEBUG: print("--------------------------------------------------------------")
                 if DEBUG: print("Fit failed for {0}. ".format(fit))
                 if DEBUG: print("--------------------------------------------------------------")
-                return None,None
+                return coeffs,1
             
                     
     

dev/incompressible_liquids/CPIncomp/DataObjects.py

@@ -207,11 +207,15 @@ class SolutionData(object):
             #if y!=0.0: raise ValueError("This is 1D only, use x not y.")
             return self.T_freeze.baseExponential(c, x)
         
+        elif self.T_freeze.type==IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL:
+            #if y!=0.0: raise ValueError("This is 1D only, use x not y.")
+            return self.T_freeze.baseLogexponential(c, x)
+        
         elif self.T_freeze.type==self.T_freeze.INCOMPRESSIBLE_EXPPOLYNOMIAL:
             return np.exp(np.polynomial.polynomial.polyval2d(p-0.0, x-self.xbase, c))
         
         else:
-            raise ValueError("Unknown function: {0}.".format(self.type))
+            raise ValueError("Unknown function: {0}.".format(self.T_freeze.type))
   
     
     

dev/incompressible_liquids/CPIncomp/ExampleObjects.py

@@ -95,7 +95,7 @@ class DigitalExamplePure(PureData,DigitalData):
         PureData.__init__(self) 
 
         self.name = "ExampleDigitalPure"
-        self.description = "water"
+        self.description = "water at 100 bar"
         self.reference = "none"
         
         self.Tmin = 280.00;

dev/incompressible_liquids/CPIncomp/SecCoolFluids.py

@@ -54,8 +54,8 @@ class SecCoolSolutionData(DigitalData):
         
         try:
             self.specific_heat.xData,self.specific_heat.yData,self.specific_heat.data = self.getArray(dataID='Cp')
-            while np.max(self.specific_heat.data[np.isfinite(self.specific_heat.data)])<500: # Expect values around 1e3
-                self.specific_heat.data *= 1e3 
+            while np.max(self.specific_heat.data[np.isfinite(self.specific_heat.data)])<1000: # Expect values around 2e3
+                self.specific_heat.data *= 10.0
             self.specific_heat.source = self.specific_heat.SOURCE_DATA
         except:
             if self.density.DEBUG: print("Could not load {}".format(self.getFile("Cp")))
@@ -63,8 +63,8 @@ class SecCoolSolutionData(DigitalData):
             
         try:
             self.conductivity.xData,self.conductivity.yData,self.conductivity.data   = self.getArray(dataID='Cond')
-            while np.max(self.conductivity.data[np.isfinite(self.conductivity.data)])>10: # Expect value below 1
-                self.conductivity.data *= 1e-3 
+            while np.max(self.conductivity.data[np.isfinite(self.conductivity.data)])>2: # Expect value below 1
+                self.conductivity.data *=  0.1 
             self.conductivity.source = self.conductivity.SOURCE_DATA
         except:
             if self.density.DEBUG: print("Could not load {}".format(self.getFile("Cond")))
@@ -72,8 +72,8 @@ class SecCoolSolutionData(DigitalData):
         
         try:
             self.viscosity.xData,self.viscosity.yData,self.viscosity.data   = self.getArray(dataID='Mu')
-            while np.max(self.viscosity.data[np.isfinite(self.viscosity.data)])>100: # Expect value below 10
-                self.viscosity.data *= 1e-3 
+            while np.max(self.viscosity.data[np.isfinite(self.viscosity.data)])>0.2: # Expect value below 0.1
+                self.viscosity.data *=  0.1
             self.viscosity.source = self.viscosity.SOURCE_DATA
         except:
             if self.density.DEBUG: print("Could not load {}".format(self.getFile("Mu")))
@@ -168,7 +168,7 @@ class SecCoolSolutionData(DigitalData):
                 if self.viscosity.coeffs==None or IncompressibleFitter.allClose(self.viscosity.coeffs, np.array([+7e+2, -6e+1, +1e+1])): # Fit failed
                     tried = True
             if len(self.viscosity.yData)>1 or tried:
-                self.viscosity.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5))
+                self.viscosity.coeffs = np.copy(std_coeffs)#np.zeros(np.round(np.array(std_coeffs.shape) * 1.5))
                 self.viscosity.type   = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL
                 self.viscosity.fitCoeffs(self.Tbase,self.xbase)
         except errList as ve:
@@ -189,13 +189,16 @@ class SecCoolSolutionData(DigitalData):
                 self.T_freeze.data = z.T
                 try:
                     self.T_freeze.source = self.T_freeze.SOURCE_DATA
-                    if np.isfinite(self.T_freeze.data).sum()<2:
-                        self.T_freeze.coeffs = np.array([+7e+6, +6e+4, +1e+1])
-                        self.T_freeze.type   = self.T_freeze.INCOMPRESSIBLE_EXPONENTIAL
-                    else:   
-                        self.T_freeze.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 2))
-                        self.T_freeze.type   = self.T_freeze.INCOMPRESSIBLE_EXPPOLYNOMIAL
+                    self.T_freeze.type   = self.T_freeze.INCOMPRESSIBLE_EXPONENTIAL
+                    self.T_freeze.coeffs = np.array([+7e+6, +6e+4, +1e+1])
                     self.T_freeze.fitCoeffs(self.Tbase,self.xbase)
+                    #if np.isfinite(self.T_freeze.data).sum()<10:
+                    #    self.T_freeze.coeffs = np.array([+7e+6, +6e+4, +1e+1])
+                    #    self.T_freeze.type   = self.T_freeze.INCOMPRESSIBLE_EXPONENTIAL
+                    #else:   
+                    #    self.T_freeze.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 2))
+                    #    self.T_freeze.type   = self.T_freeze.INCOMPRESSIBLE_EXPPOLYNOMIAL
+                    #self.T_freeze.fitCoeffs(self.Tbase,self.xbase)
                 except errList as ve:
                     if self.T_freeze.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"T_freeze",ve))
                     pass

dev/incompressible_liquids/CPIncomp/WriterObjects.py

@@ -8,6 +8,7 @@ from CPIncomp.DataObjects import SolutionData
 from CPIncomp.BaseObjects import IncompressibleData, IncompressibleFitter
 from matplotlib.patches import Rectangle
 from matplotlib.ticker import MaxNLocator
+from matplotlib.backends.backend_pdf import PdfPages
 
 class SolutionDataWriter(object):
     """ 
@@ -76,7 +77,8 @@ class SolutionDataWriter(object):
                 if fluidObject.viscosity.coeffs==None or IncompressibleFitter.allClose(fluidObject.viscosity.coeffs, np.array([+5e+2, -6e+1, +1e+1])): # Fit failed
                     tried = True
             if len(fluidObject.viscosity.yData)>1 or tried:
-                fluidObject.viscosity.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5))
+                #fluidObject.viscosity.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5))
+                fluidObject.viscosity.coeffs = np.copy(std_coeffs)
                 fluidObject.viscosity.type   = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL
                 fluidObject.viscosity.fitCoeffs(tBase,xBase)
         except errList as ve:
@@ -93,7 +95,8 @@ class SolutionDataWriter(object):
                 if fluidObject.saturation_pressure.coeffs==None or IncompressibleFitter.allClose(fluidObject.saturation_pressure.coeffs, np.array([-5e+3, +6e+1, -1e+1])): # Fit failed
                     tried = True
             if len(fluidObject.saturation_pressure.yData)>1 or tried:
-                fluidObject.saturation_pressure.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5))
+                #fluidObject.saturation_pressure.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5))
+                fluidObject.saturation_pressure.coeffs = np.copy(std_coeffs)
                 fluidObject.saturation_pressure.type   = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL
                 fluidObject.saturation_pressure.fitCoeffs(tBase,xBase)
         except errList as ve:
@@ -104,11 +107,11 @@ class SolutionDataWriter(object):
         if fluidObject.xid!=fluidObject.ifrac_pure and fluidObject.xid!=fluidObject.ifrac_undefined:
             fluidObject.T_freeze.setxyData([0.0],xData)
             try:
-                if len(fluidObject.T_freeze.xData)==1 and np.isfinite(fluidObject.T_freeze.data).sum()<2:
+                if len(fluidObject.T_freeze.xData)==1:# and np.isfinite(fluidObject.T_freeze.data).sum()<10:
                     fluidObject.T_freeze.coeffs = np.array([+7e+2, -6e+1, +1e+1])
                     fluidObject.T_freeze.type   = IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL
                 else:   
-                    fluidObject.T_freeze.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1))
+                    fluidObject.specific_heat.coeffs = np.copy(std_coeffs)
                     fluidObject.T_freeze.type   = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL
                 fluidObject.T_freeze.fitCoeffs(tBase,xBase)
             except errList as ve:
@@ -276,11 +279,13 @@ class SolutionDataWriter(object):
         print(" ... done") 
         return 
     
-    def writeReportList(self, fluidObjs):
+    def writeReportList(self, fluidObjs, pdfFile=None):
         print("Writing fitting reports:", end="")
+        pdfObj = None
+        if pdfFile!=None: pdfObj = PdfPages(pdfFile)
         for obj in fluidObjs:
             self.printStatusID(fluidObjs, obj)
-            self.makeFitReportPage(obj)           
+            self.makeFitReportPage(obj,pdfObj=pdfObj)           
             try: 
                 self.makeFitReportPage(obj)
             except (TypeError, ValueError) as e:
@@ -288,6 +293,7 @@ class SolutionDataWriter(object):
                 print(obj)
                 print(e)
                 pass
+        if pdfFile!=None: pdfObj.close()
         print(" ... done") 
         return 
     
@@ -552,6 +558,10 @@ class SolutionDataWriter(object):
 
         if zFunc!=None and axVal!=None:
             axVal.plot(pFunc, zFunc, label='function' , **fitFormatter)
+            if solObj.xid!=solObj.ifrac_pure and not xFunction:
+                axVal.set_title("showing x={0:3.2f}".format(xFunc[0]))
+            else:
+                axVal.set_title(" ")
             
         if zMiMa!=None and axVal!=None:
             axVal.plot(pMiMa, zMiMa, alpha=0.25, ls=':', color=fitFormatter["color"])
@@ -561,10 +571,7 @@ class SolutionDataWriter(object):
             
         if zError!=None and axErr!=None:
             axErr.plot(pData, zError, label='error' , **errorFormatter)
-            if solObj.xid!=solObj.ifrac_pure and not xFunction:
-                axErr.set_title("showing x={0:3.2f}".format(xData[0]))
-            else:
-                axErr.set_title(" ")
+
         elif axErr!=None:
             errorFormatter['alpha'] = 0.00
             axErr.plot([pData[0],pData[-1]], [0,0], **errorFormatter)
@@ -642,29 +649,42 @@ class SolutionDataWriter(object):
         if solObj.xid==solObj.ifrac_volume: conc=True
         if solObj.xid==solObj.ifrac_mole: conc=True
         if conc==True:
-            myAnnotate('Composition: ',u'{0} % to {1} %, {2}'.format(solObj.xmin*100., solObj.xmax*100., solObj.xid),x=x,y=y); x += .0; y -= dy
+            myAnnotate('Composition: ',u'{0} % to {1} %, {2}'.format(solObj.xmin*100., solObj.xmax*100., solObj.xid),x=x,y=y)
         else:
-            x += .0; y -= dy
+            myAnnotate('Composition: ','pure fluid',x=x,y=y)
+        x += .0; y -= dy
             
         if solObj.density.source!=solObj.density.SOURCE_NOT_SET:
             myAnnotate('Density: ',u'{0} to {1} {2}'.format(solObj.density.source, solObj.density.type, solObj.density.coeffs.shape),x=x,y=y)
+        else:
+            myAnnotate('Density: ','no information',x=x,y=y)
         x += .0; y -= dy
         if solObj.specific_heat.source!=solObj.specific_heat.SOURCE_NOT_SET:
             myAnnotate('Spec. Heat: ',u'{0} to {1} {2}'.format(solObj.specific_heat.source, solObj.specific_heat.type, solObj.specific_heat.coeffs.shape),x=x,y=y)
+        else:
+            myAnnotate('Spec. Heat: ','no information',x=x,y=y)
         x += .0; y -= dy
             
         x = xStart + dx; y = yStart-dy-dy
         if solObj.conductivity.source!=solObj.conductivity.SOURCE_NOT_SET:
             myAnnotate('Th. Cond.: ',u'{0} to {1} {2}'.format(solObj.conductivity.source, solObj.conductivity.type, solObj.conductivity.coeffs.shape),x=x,y=y)
+        else:
+            myAnnotate('Th. Cond.: ','no information',x=x,y=y)
         x += .0; y -= dy
         if solObj.viscosity.source!=solObj.viscosity.SOURCE_NOT_SET:
             myAnnotate('Viscosity: ',u'{0} to {1} {2}'.format(solObj.viscosity.source, solObj.viscosity.type, solObj.viscosity.coeffs.shape),x=x,y=y)
+        else:
+            myAnnotate('Viscosity: ','no information',x=x,y=y)
         x += .0; y -= dy
         if solObj.saturation_pressure.source!=solObj.saturation_pressure.SOURCE_NOT_SET:
             myAnnotate('Psat: ',u'{0} to {1} {2}'.format(solObj.saturation_pressure.source, solObj.saturation_pressure.type, solObj.saturation_pressure.coeffs.shape),x=x,y=y)
+        else:
+            myAnnotate('Psat: ','no information',x=x,y=y)
         x += .0; y -= dy
         if solObj.T_freeze.source!=solObj.T_freeze.SOURCE_NOT_SET:
             myAnnotate('Tfreeze: ',u'{0} to {1} {2}'.format(solObj.T_freeze.source, solObj.T_freeze.type, solObj.T_freeze.coeffs.shape),x=x,y=y)
+        else:
+            myAnnotate('Tfreeze: ','no information',x=x,y=y)
         x += .0; y -= dy
 
         
@@ -682,7 +702,7 @@ class SolutionDataWriter(object):
     
     
     
-    def makeFitReportPage(self, solObj=SolutionData()):
+    def makeFitReportPage(self, solObj=SolutionData(), pdfObj=None):
         """
         Creates a whole page with some plots and basic information
         for both fit quality, reference data, data sources and 
@@ -866,11 +886,8 @@ class SolutionDataWriter(object):
         #plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
         
         plt.savefig(os.path.join("report","{0}_fitreport.pdf".format(solObj.name)))
+        if pdfObj!=None: pdfObj.savefig(fig)
         plt.close()
-
-        
-        
-        
         pass
     
     

dev/incompressible_liquids/all_incompressibles.py

@@ -3,6 +3,7 @@ from CPIncomp.WriterObjects import SolutionDataWriter
 from CPIncomp import getExampleNames, getPureFluids, getCoefficientFluids,\
     getDigitalFluids, getSecCoolFluids, getMelinderFluids
 import sys
+from CPIncomp.DataObjects import SolutionData
 
 
 if __name__ == '__main__':   
@@ -10,28 +11,28 @@ if __name__ == '__main__':
     writer = SolutionDataWriter()
     doneObjs = []
     
-    # To debug single fluids
-    from CPIncomp.SecCoolFluids import SecCoolSolutionData,SecCoolIceData
-    from CPIncomp.PureFluids import Therminol72
-    #solObjs  = [SecCoolSolutionData(sFile='Melinder, Ammonia'            ,sFolder='xMass',name='MAM2',desc='Melinder, Ammonia'            ,ref='Melinder-BOOK-2010, SecCool software')]
-    #solObjs += [SecCoolIceData(sFile='IceNA'   ,sFolder='xMass',name='IceNA',desc='Ice slurry with NaCl' ,ref='Danish Technological Institute, SecCool software')]
-    #solObjs = [Freezium()]
-    #solObjs[0].density.DEBUG = True
-    #solObjs[0].specific_heat.DEBUG = True
-    #solObjs[0].conductivity.DEBUG = True
-    #solObjs[0].viscosity.DEBUG = True
-    #solObjs[0].T_freeze.DEBUG = True
-    #writer.fitSecCoolList(solObjs)
-    solObjs = [Therminol72()]
-    solObjs[0].viscosity.DEBUG=True
-    #solObjs[0].saturation_pressure.DEBUG=True
-    writer.fitFluidList([solObjs[-1]])
-    #
-    ##from CPIncomp.ExampleObjects import SecCoolExample
-    ##solObjs = [SecCoolExample()]
-    writer.writeFluidList(solObjs)
-    writer.writeReportList(solObjs)
-    sys.exit(0)
+#    # To debug single fluids
+#    from CPIncomp.SecCoolFluids import SecCoolSolutionData,SecCoolIceData
+#    from CPIncomp.PureFluids import Therminol72, DowthermJ
+#    #solObjs  = [SecCoolSolutionData(sFile='Melinder, Ammonia'            ,sFolder='xMass',name='MAM2',desc='Melinder, Ammonia'            ,ref='Melinder-BOOK-2010, SecCool software')]
+#    #solObjs += [SecCoolIceData(sFile='IceNA'   ,sFolder='xMass',name='IceNA',desc='Ice slurry with NaCl' ,ref='Danish Technological Institute, SecCool software')]
+#    #solObjs = [Freezium()]
+#    #solObjs[0].density.DEBUG = True
+#    #solObjs[0].specific_heat.DEBUG = True
+#    #solObjs[0].conductivity.DEBUG = True
+#    #solObjs[0].viscosity.DEBUG = True
+#    #solObjs[0].T_freeze.DEBUG = True
+#    #writer.fitSecCoolList(solObjs)
+#    solObjs = [Therminol72()]#,Therminol72()]
+#    #solObjs[0].viscosity.DEBUG=True
+#    solObjs[0].saturation_pressure.DEBUG=True
+#    #
+#    ##from CPIncomp.ExampleObjects import SecCoolExample
+#    ##solObjs = [SecCoolExample()]
+#    writer.fitFluidList(solObjs)
+#    writer.writeFluidList(solObjs)
+#    writer.writeReportList(solObjs)
+#    sys.exit(0)
     
     fluidObjs = getExampleNames(obj=True)
     examplesToFit = ["ExamplePure","ExampleSolution","ExampleDigital","ExampleDigitalPure"]
@@ -43,7 +44,7 @@ if __name__ == '__main__':
     print("\nProcessing example fluids")
     writer.writeFluidList(doneObjs)
     
-    writer.writeReportList(doneObjs)
+    writer.writeReportList(doneObjs, pdfFile="all_examples.pdf")
     #sys.exit(0)
         
     # If the examples did not cause any errors, 
@@ -77,16 +78,39 @@ if __name__ == '__main__':
     print("Checking the list of fluid objects.")
     #doneObjs += getCoefficientObjects()[:]
     doneObjs = sorted(doneObjs, key=lambda x: x.name)
+    
+    purefluids= []
+    solutions = []
+    errors    = []
     for i in range(len(doneObjs)-1):
         if doneObjs[i].name==doneObjs[i+1].name:
             print("Conflict between {0} and {1}, aborting".format(doneObjs[i],doneObjs[i+1]))
             raise ValueError("Two elements have the same name, that does not work: {0}".format(doneObjs[i].name))
+        else:
+            if doneObjs[i].xid==SolutionData.ifrac_mass or \
+              doneObjs[i].xid==SolutionData.ifrac_mole or \
+              doneObjs[i].xid==SolutionData.ifrac_volume:
+                solutions += [doneObjs[i]]
+            elif doneObjs[i].xid==SolutionData.ifrac_pure:
+                purefluids += [doneObjs[i]]
+            else: errors += [doneObjs[i]]
+        
+    print("All checks passed, going to write parameters to disk.")
+    writer.writeFluidList(doneObjs)    
     
-    
-    print("All checks passed, going to write to disk.")
-    
-    writer.writeFluidList(doneObjs)
-    writer.writeReportList(doneObjs)
+    print("Creating the fitting reports for the different groups.")
+    #writer.writeReportList(doneObjs)
+    #doneObjs.sort(key=lambda x: (x.xid ,x.name))
+    if len(purefluids)>0:
+        print("Processing {0:2d} pure fluids   - ".format(len(purefluids)), end="")
+        writer.writeReportList(purefluids, pdfFile="all_pure.pdf")
+    if len(solutions)>0:
+        print("Processing {0:2d} solutions     - ".format(len(solutions)), end="")
+        writer.writeReportList(solutions, pdfFile="all_solutions.pdf")  
+    if len(errors)>0:
+        print("Processing {0:2d} faulty fluids - ".format(len(errors)), end="")
+        writer.writeReportList(errors, pdfFile="all_errors.pdf")
+       
 
     print("All done, bye")
     sys.exit(0)

include/DataStructures.h

@@ -74,10 +74,10 @@ int get_parameter_index(const std::string &param_name);
 std::string get_csv_parameter_list();
 
 /// These are constants for the compositions
-enum compositions {ifrac_mass, ifrac_mole, ifrac_volume, ifrac_undefined, ifrac_pure};
+enum compositions{ifrac_mass, ifrac_mole, ifrac_volume, ifrac_undefined, ifrac_pure};
 
 /// These are constants for the phases of the fluid
-enum phases {iphase_liquid, iphase_supercritical, iphase_gas, iphase_twophase, iphase_unknown};
+enum phases{ iphase_liquid, iphase_supercritical, iphase_gas, iphase_twophase, iphase_unknown };
 
 /// These are unit types for the fluid
 enum fluid_types{FLUID_TYPE_PURE, FLUID_TYPE_PSEUDOPURE, FLUID_TYPE_REFPROP, FLUID_TYPE_INCOMPRESSIBLE_LIQUID, FLUID_TYPE_INCOMPRESSIBLE_SOLUTION, FLUID_TYPE_UNDEFINED};

include/IncompressibleFluid.h

@@ -32,6 +32,7 @@ struct IncompressibleData {
 		INCOMPRESSIBLE_POLYNOMIAL,
 		INCOMPRESSIBLE_EXPPOLYNOMIAL,
 		INCOMPRESSIBLE_EXPONENTIAL,
+		INCOMPRESSIBLE_LOGEXPONENTIAL,
 		INCOMPRESSIBLE_POLYOFFSET
 	};
 	Eigen::MatrixXd coeffs; //TODO: Can we store the Eigen::Matrix objects more efficiently?
@@ -191,6 +192,7 @@ public:
 protected:
 	/// Base functions that handle the custom function types
 	double baseExponential(IncompressibleData data, double y, double ybase);
+	double baseLogexponential(IncompressibleData data, double y, double ybase);
 	double baseExponentialOffset(IncompressibleData data, double y);
 	double basePolyOffset(IncompressibleData data, double y, double z=0.0);
 
@@ -309,6 +311,7 @@ protected:
 	 *  xmax since the default values cause an error. */
 	bool checkX(double x);
 
+public:
 	/// Check validity of temperature, pressure and composition input.
 	bool checkTPX(double T, double p, double x){
 		return (checkT(T,p,x) && checkP(T,p,x) && checkX(x));

src/Backends/Incompressible/IncompressibleBackend.cpp

@@ -101,6 +101,7 @@ void IncompressibleBackend::update(long input_pair, double value1, double value2
     if (_T < 0){ throw ValueError("T is less than zero");}
     if (!ValidNumber(_T)){ throw ValueError("T is not a valid number");}
     if (get_debug_level()>=50) std::cout << format("Incompressible backend: Update finished T=%f, p=%f, x=%s ",this->_T,this->_p,vec_to_string(_fractions).c_str()) << std::endl;
+    fluid->checkTPX(_T,_p,_fractions[0]);
 }
 
 /// Set the mole fractions
@@ -468,6 +469,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
 		CAPTURE(res);
 		CHECK( check_abs(val,res,acc) );
 		}
+		CoolProp::set_debug_level(100);
 		// ... as %
 		res = CoolProp::PropsSI("D","T",T,"P",p,fluid+format("-%f%s",x*100.0,"%"));
 		{

src/Backends/Incompressible/IncompressibleBackend.h

@@ -2,9 +2,9 @@
 #ifndef INCOMPRESSIBLEBACKEND_H_
 #define INCOMPRESSIBLEBACKEND_H_
 
+#include "DataStructures.h"
 #include "IncompressibleFluid.h"
 #include "AbstractState.h"
-#include "DataStructures.h"
 #include "Exceptions.h"
 
 #include <vector>
@@ -33,9 +33,9 @@ public:
 	IncompressibleBackend(const std::vector<std::string> &component_names);
 
     // Incompressible backend uses different compositions
-    bool using_mole_fractions(){return this->fluid->getxid()==ifrac_mole;};
+    bool using_mole_fractions(){return  this->fluid->getxid()==ifrac_mole;};
     bool using_mass_fractions(){return (this->fluid->getxid()==ifrac_mass || this->fluid->getxid()==ifrac_pure);};
-    bool using_volu_fractions(){return this->fluid->getxid()==ifrac_volume;};
+    bool using_volu_fractions(){return  this->fluid->getxid()==ifrac_volume;};
 
     /// Updating function for incompressible fluid
     /**

src/Backends/Incompressible/IncompressibleFluid.cpp

@@ -1,9 +1,9 @@
 
+#include "DataStructures.h"
 #include "IncompressibleFluid.h"
 #include "math.h"
 #include "MatrixMath.h"
 #include "PolyMath.h"
-#include "DataStructures.h"
 
 namespace CoolProp {
 
@@ -39,13 +39,21 @@ bool IncompressibleFluid::is_pure() {
 	return false;
 }
 
-/// Base functions that handle the custom data type, just a place holder to show the structure.
+/// Base exponential function
 double IncompressibleFluid::baseExponential(IncompressibleData data, double y, double ybase){
 	Eigen::VectorXd coeffs = makeVector(data.coeffs);
 	size_t r=coeffs.rows(),c=coeffs.cols();
 	if (strict && (r!=3 || c!=1) ) throw ValueError(format("%s (%d): You have to provide a 3,1 matrix of coefficients, not  (%d,%d).",__FILE__,__LINE__,r,c));
 	return exp( (double) (coeffs[0] / ( (y-ybase)+coeffs[1] ) - coeffs[2] ) );
 }
+/// Base exponential function with logarithmic term
+double IncompressibleFluid::baseLogexponential(IncompressibleData data, double y, double ybase){
+	Eigen::VectorXd coeffs = makeVector(data.coeffs);
+	size_t r=coeffs.rows(),c=coeffs.cols();
+	if (strict && (r!=3 || c!=1) ) throw ValueError(format("%s (%d): You have to provide a 3,1 matrix of coefficients, not  (%d,%d).",__FILE__,__LINE__,r,c));
+	return exp( (double) ( log( (double) (1.0/((y-ybase)+coeffs[0]) + 1.0/((y-ybase)+coeffs[0])/((y-ybase)+coeffs[0]) ) ) *coeffs[1]+coeffs[2] ) );
+}
+
 double IncompressibleFluid::basePolyOffset(IncompressibleData data, double y, double z){
 	size_t r=data.coeffs.rows(),c=data.coeffs.cols();
 	double offset = 0.0;
@@ -77,6 +85,9 @@ double IncompressibleFluid::rho (double T, double p, double x){
 		case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
 			return baseExponential(density, T, 0.0);
 			break;
+		case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+			return baseLogexponential(density, T, 0.0);
+			break;
 		case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 			return exp(poly.evaluate(density.coeffs, T, x, 0, 0, Tbase, xbase));
 			break;
@@ -156,6 +167,9 @@ double IncompressibleFluid::visc(double T, double p, double x){
 		case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
 			return baseExponential(viscosity, T, 0.0);
 			break;
+		case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+			return baseLogexponential(viscosity, T, 0.0);
+			break;
 		case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 			return exp(poly.evaluate(viscosity.coeffs, T, x, 0, 0, Tbase, xbase));
 			break;
@@ -180,6 +194,9 @@ double IncompressibleFluid::cond(double T, double p, double x){
 		case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
 			return baseExponential(conductivity, T, 0.0);
 			break;
+		case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+			return baseLogexponential(conductivity, T, 0.0);
+			break;
 		case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 			return exp(poly.evaluate(conductivity.coeffs, T, x, 0, 0, Tbase, xbase));
 			break;
@@ -205,6 +222,9 @@ double IncompressibleFluid::psat(double T,           double x){
 		case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
 			return baseExponential(p_sat, T, 0.0);
 			break;
+		case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+			return baseLogexponential(p_sat, T, 0.0);
+			break;
 		case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 			return exp(poly.evaluate(p_sat.coeffs, T, x, 0, 0, Tbase, xbase));
 			break;
@@ -229,6 +249,9 @@ double IncompressibleFluid::Tfreeze(       double p, double x){
 		case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
 			return baseExponential(T_freeze, x, 0.0);
 			break;
+		case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+			return baseLogexponential(T_freeze, x, 0.0);
+			break;
 		case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 			return exp(poly.evaluate(T_freeze.coeffs, p, x, 0, 0, 0.0, xbase));
 			break;
@@ -260,7 +283,10 @@ double IncompressibleFluid::inputFromMass (double T,     double x){
 				return poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way
 				break;
 			case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
-				return baseExponential(mass2input, T, 0.0);
+				return baseExponential(mass2input, x, 0.0);
+				break;
+			case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+				return baseLogexponential(mass2input, x, 0.0);
 				break;
 			case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 				return exp(poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way
@@ -294,7 +320,10 @@ double IncompressibleFluid::inputFromVolume (double T,   double x){
 				return poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way
 				break;
 			case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
-				return baseExponential(volume2input, T, 0.0);
+				return baseExponential(volume2input, x, 0.0);
+				break;
+			case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+				return baseLogexponential(volume2input, x, 0.0);
 				break;
 			case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 				return exp(poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way
@@ -328,7 +357,10 @@ double IncompressibleFluid::inputFromMole (double T,     double x){
 				return poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way
 				break;
 			case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
-				return baseExponential(mole2input, T, 0.0);
+				return baseExponential(mole2input, x, 0.0);
+				break;
+			case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
+				return baseLogexponential(mole2input, x, 0.0);
 				break;
 			case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
 				return exp(poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way
@@ -441,21 +473,11 @@ double IncompressibleFluid::T_u   (double Umass, double p, double x){
 bool IncompressibleFluid::checkT(double T, double p, double x) {
 	if (Tmin <= 0.) throw ValueError("Please specify the minimum temperature.");
 	if (Tmax <= 0.) throw ValueError("Please specify the maximum temperature.");
-	if ((Tmin > T) || (T > Tmax)) throw ValueError(
-				format("Your temperature %f is not between %f and %f.",
-						T, Tmin, Tmax));
+	if ((Tmin > T) || (T > Tmax)) throw ValueError(format("Your temperature %f is not between %f and %f.", T, Tmin, Tmax));
 	double TF = 0.0;
-	if (T_freeze.type!=IncompressibleData::INCOMPRESSIBLE_NOT_SET) {
-		TF = Tfreeze(p, x);
-	}
-	if ( T<TF) {
-		throw ValueError(
-			format("Your temperature %f is below the freezing point of %f.",
-					T, TF));
-	} else {
-		return true;
-	}
-	return false;
+	if (T_freeze.type!=IncompressibleData::INCOMPRESSIBLE_NOT_SET) TF = Tfreeze(p, x);
+	if ( T<TF) throw ValueError(format("Your temperature %f is below the freezing point of %f.", T, TF));
+	return true;
 }
 
 /// Check validity of pressure input.
@@ -467,21 +489,10 @@ bool IncompressibleFluid::checkT(double T, double p, double x) {
  *  */
 bool IncompressibleFluid::checkP(double T, double p, double x) {
 	double ps = 0.0;
-	if (p_sat.type!=IncompressibleData::INCOMPRESSIBLE_NOT_SET) {
-		ps = psat(T, x);
-	}
-	if (p < 0.0) {
-		throw ValueError(
-				format("You cannot use negative pressures: %f < %f. ",
-						p, 0.0));
-	} else if (p < ps) {
-		throw ValueError(
-				format("Equations are valid for liquid phase only: %f < %f. ",
-						p, ps));
-	} else {
-		return true;
-	}
-	return false;
+	if (p_sat.type!=IncompressibleData::INCOMPRESSIBLE_NOT_SET) ps = psat(T, x);
+	if (p < 0.0) throw ValueError(format("You cannot use negative pressures: %f < %f. ", p, 0.0));
+	if (ps>0.0 && p < ps)  throw ValueError(format("Equations are valid for liquid phase only: %f < %f (psat). ", p, ps));
+	return true;
 }
 
 /// Check validity of composition input.
@@ -489,18 +500,10 @@ bool IncompressibleFluid::checkP(double T, double p, double x) {
  *  maximum value. Enforces the redefinition of xmin and
  *  xmax since the default values cause an error. */
 bool IncompressibleFluid::checkX(double x){
-	if (xmin < 0. || xmin > 1.0) {
-		throw ValueError("Please specify the minimum concentration between 0 and 1.");
-	} else if (xmax < 0.0 || xmax > 1.0) {
-		throw ValueError("Please specify the maximum concentration between 0 and 1.");
-	} else if ((xmin > x) || (x > xmax)) {
-		throw ValueError(
-				format("Your composition %f is not between %f and %f.",
-						x, xmin, xmax));
-	} else {
-		return true;
-	}
-	return false;
+	if (xmin < 0.0 || xmin > 1.0) throw ValueError("Please specify the minimum concentration between 0 and 1.");
+	if (xmax < 0.0 || xmax > 1.0) throw ValueError("Please specify the maximum concentration between 0 and 1.");
+	if ((xmin > x) || (x > xmax)) throw ValueError(format("Your composition %f is not between %f and %f.", x, xmin, xmax));
+	return true;
 }
 
 } /* namespace CoolProp */

src/Backends/Incompressible/IncompressibleLibrary.cpp

@@ -367,6 +367,11 @@ IncompressibleData JSONIncompressibleLibrary::parse_coefficients(rapidjson::Valu
 					fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
 					return fluidData;
 				}
+				else if (!type.compare("logexponential")){
+					fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL;
+					fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
+					return fluidData;
+				}
 				else if (!type.compare("exppolynomial")){
 					fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL;
 					fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));