fits, fits, fits

dev/incompressible_liquids/CPIncomp/BaseObjects.py

@@ -17,6 +17,7 @@ class IncompressibleData(object):
     INCOMPRESSIBLE_POLYNOMIAL    = 'polynomial'
     INCOMPRESSIBLE_EXPPOLYNOMIAL = 'exppolynomial'
     INCOMPRESSIBLE_EXPONENTIAL   = 'exponential'
+    INCOMPRESSIBLE_LOGEXPONENTIAL= 'logexponential'
     INCOMPRESSIBLE_POLYOFFSET    = 'polyoffset'
     INCOMPRESSIBLE_CHEBYSHEV     = 'chebyshev'
     
@@ -35,6 +36,7 @@ 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
         
     @staticmethod
@@ -53,6 +55,10 @@ class IncompressibleData(object):
         elif eqnType==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL:
             #if y!=0.0: raise ValueError("This is 1D only, use x not y.")
             return IncompressibleData.baseExponential(c, x)
+            
+        elif eqnType==IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL:
+            #if y!=0.0: raise ValueError("This is 1D only, use x not y.")
+            return IncompressibleData.baseLogexponential(c, x)
         
         elif eqnType==IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL:
             return np.exp(np.polynomial.polynomial.polyval2d(x-xbase, y-ybase, c))
@@ -67,6 +73,14 @@ class IncompressibleData(object):
             raise ValueError("You have to provide a 3,1 matrix of coefficients, not ({0},{1}).".format(r,c))
         coeffs_tmp = np.array(coeffs.flat)
         return np.exp(np.clip( (coeffs_tmp[0]/ ( x+coeffs_tmp[1] ) - coeffs_tmp[2]),IncompressibleData.minLog,IncompressibleData.maxLog))
+        
+    @staticmethod
+    def baseLogexponential(co, x):
+        r,c,coeffs = IncompressibleFitter.shapeArray(co)
+        if not ( (r==4 and c==1) or (r==1 and c==4) ):
+            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))
 
     @staticmethod
     def basePolyOffset(co, x):
@@ -85,25 +99,51 @@ class IncompressibleData(object):
     
     def fitCoeffs(self, xbase, ybase, x=None, y=None):
         
-        if (x!=None and self.xData!=None and not np.allclose(x, self.xData)) \
+        if (x!=None and self.xData!=None and not IncompressibleFitter.allClose(x, self.xData)) \
         or (x==None and self.xData==None): raise ValueError("I do not know which x-value you would like to use. Define either x or self.xData.")
-        if (y!=None and self.yData!=None and not np.allclose(y, self.yData)) \
+        if (y!=None and self.yData!=None and not IncompressibleFitter.allClose(y, self.yData)) \
         or (y==None and self.yData==None): raise ValueError("I do not know which y-value you would like to use. Define either y or self.yData.")
         
         if x==None and self.xData!=None: x=self.xData
         if y==None and self.yData!=None: y=self.yData
         
-        res = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \
+        #res = None
+        #r2 = None
+        
+        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 = 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)
+            
+            #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
+         
         if res==None:
-            raise ValueError("There was a fitting error.")
-        elif np.allclose(res, self.coeffs):
+            if self.DEBUG: print("There was a fitting error, no solution found.")
+        elif IncompressibleFitter.allClose(res, self.coeffs):
             if self.DEBUG: print("Coefficients did not change.")
         else:
             self.coeffs = res
+            self.rsq    = r2 
             
 
     def setxData(self, xData):
@@ -138,6 +178,12 @@ class IncompressibleData(object):
 
 class IncompressibleFitter(object):
     
+    @staticmethod
+    def allClose(a,b):
+        if np.array(a).shape==np.array(b).shape:
+            return np.allclose(a, b)
+        else: return False
+    
     @staticmethod
     def shapeArray(array, axs=0):
         """
@@ -239,7 +285,7 @@ class IncompressibleFitter(object):
             if cr==1 and cc==1: 
                 if DEBUG: print("No coefficients left to fit, aborting procedure.")
                 coeffs = np.array([[0.0]])
-                return
+                return coeffs, None 
             if (xr<cr): 
                 if DEBUG: print("Less data points than coefficients in first dimension ({0} < {1}), reducing coefficient matrix.".format(xr,cr))
                 coeffs = coeffs[:xr,:]
@@ -252,9 +298,10 @@ class IncompressibleFitter(object):
             z_input = np.copy(z)
             if eqnType==IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL:
                 z_input = np.log(z_input)
-            coeffs = IncompressibleFitter.getCoeffs2d(x_input, y_input, z_input, cr-1, cc-1, DEBUG=DEBUG)
+            coeffs,r2 = IncompressibleFitter.getCoeffs2d(x_input, y_input, z_input, cr-1, cc-1, DEBUG=DEBUG)
             if DEBUG: print("Coefficients after fitting: \n{0}".format(coeffs))
-            return coeffs
+            if DEBUG: print("Coefficient of determination: {0}".format(r2))
+            return coeffs,r2
         
         
         # Select if 1D or 2D fitting
@@ -262,13 +309,14 @@ class IncompressibleFitter(object):
             if DEBUG: print("1D function detected.")
             if yc==1: 
                 if DEBUG: print("Fitting {0} in x-direction.".format(eqnType))
-                coeffs = IncompressibleFitter.getCoeffsIterative1D(x, z, eqnType=eqnType, coeffs=coeffs, DEBUG=DEBUG)
+                coeffs,r2 = IncompressibleFitter.getCoeffsIterative1D(x, z, eqnType=eqnType, coeffs=coeffs, DEBUG=DEBUG)
             elif xr==1: 
                 if DEBUG: print("Fitting {0} in y-direction.".format(eqnType))
-                coeffs = IncompressibleFitter.getCoeffsIterative1D(y, z, eqnType=eqnType, coeffs=coeffs, DEBUG=DEBUG)
+                coeffs,r2 = 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))
-            return coeffs
+            if DEBUG: print("Coefficient of determination: {0}".format(r2))
+            return coeffs, r2
 
         elif yc>1: # 2D fitting
             raise ValueError("There are no other 2D fitting functions than polynomials, cannot use {0}.".format(eqnType))
@@ -297,6 +345,8 @@ 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
@@ -349,12 +399,18 @@ 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
+        #print("\n r2 2d: ",r2.shape,r2,"\n")
+        
         #Rearrange coefficients to a matrix shape
         C = np.zeros((x_order,y_order))
         for i, (xi,yi) in enumerate(xy_exp): # makes columns
             C[xi][yi] = coeffs[i]
         
-        return C
+        return C, r2
     
     @staticmethod
     def getCoeffsIterative1D(x_in, z_in, eqnType, coeffs, DEBUG=False):
@@ -365,6 +421,8 @@ 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
+        
         #fit = "Powell" # use Powell's algorithm
         #fit = "BFGS" # use Broyden-Fletcher-Goldfarb-Shanno
         #fit = "LMA" # use the Levenberg-Marquardt algorithm from curve_fit 
@@ -377,7 +435,7 @@ class IncompressibleFitter(object):
         
         expLog = False
         # Fitting the logarithm of z_in?
-        if (eqnType==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL):
+        if (eqnType==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL or eqnType==IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL):
             expLog = True
         
         xData = np.array(x_in.flat)
@@ -429,16 +487,20 @@ class IncompressibleFitter(object):
 #                        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: ".format(pcov))
-                        return 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 )
+                        #print("\n r2 FMA: ",r2.shape,r2,"\n")
+                        return popt,r2
                     else:
-                        print("Fit failed for {0}.".format(algorithm))
-                        sys.stdout.flush()
+                        if DEBUG: print("Fit failed for {0}.".format(algorithm))
+                        if DEBUG: sys.stdout.flush()
                         success = False 
                     
                 except RuntimeError as e:
-                    print("Exception using "+algorithm+": "+str(e))
-                    sys.stdout.flush()
+                    if DEBUG: print("Exception using "+algorithm+": "+str(e))
+                    if DEBUG: sys.stdout.flush()
                     success = False
         
             #fit = "MIN" # use a home-made minimisation with Powell and Broyden-Fletcher-Goldfarb-Shanno
@@ -452,14 +514,19 @@ class IncompressibleFitter(object):
                     if res.success:
                         success = True
                         if DEBUG: print("Fit succeeded with: {0}".format(algorithm))
-                        return res.x
+                        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
                     else:
-                        print("Fit failed for {0}.".format(algorithm))
-                        sys.stdout.flush()
+                        if DEBUG: print("Fit failed for {0}.".format(algorithm))
+                        if DEBUG: sys.stdout.flush()
                         success = False
                 except RuntimeError as e:
-                    print("Exception using "+algorithm+": "+str(e))
-                    sys.stdout.flush()
+                    if DEBUG: print("Exception using "+algorithm+": "+str(e))
+                    if DEBUG: sys.stdout.flush()
                     success = False
 
             # Something went wrong, probably a typo in the algorithm selector
@@ -470,16 +537,16 @@ class IncompressibleFitter(object):
                 #print("Fit did not succeed with {0}, reducing tolerance to {1}.".format(algorithm,tol))
                 success = False
                 counter += 1
-            elif tolerance<1e-8:
+            elif tolerance<1e-10:
                 tolerance *= 1e2
-                print("Fit did not succeed, reducing tolerance to {0}.".format(tolerance))
+                if DEBUG: print("Fit did not succeed, reducing tolerance to {0}.".format(tolerance))
                 success = False
                 counter = 0                
             else:
-                print("--------------------------------------------------------------")
-                print("Fit failed for {0}. ".format(fit))
-                print("--------------------------------------------------------------")
-                return None 
+                if DEBUG: print("--------------------------------------------------------------")
+                if DEBUG: print("Fit failed for {0}. ".format(fit))
+                if DEBUG: print("--------------------------------------------------------------")
+                return None,None
             
                     
     

dev/incompressible_liquids/CPIncomp/DataObjects.py

@@ -108,9 +108,45 @@ class SolutionData(object):
 #        objList["viscosity"] = self.viscosity
 #        objList["saturation pressure"] = self.saturation_pressure
 #        return objList
+
+
+    def checkT(self, T, p, x):
+        if self.Tmin <= 0.: raise ValueError("Please specify the minimum temperature.")
+        if self.Tmax <= 0.: raise ValueError("Please specify the maximum temperature.");
+        if ((self.Tmin > T) or (T > self.Tmax)): raise ValueError("Your temperature {0} is not between {1} and {2}.".format(T, self.Tmin, self.Tmax))
+        TF = 0.0
+        if (self.T_freeze.type!=IncompressibleData.INCOMPRESSIBLE_NOT_SET): TF = self.Tfreeze(T, p, x)
+        if ( T<TF ): raise ValueError("Your temperature {0} is below the freezing point of {1}.".format(T, TF))
+        else: return True
+        return False
+
+    def checkP(self, T, p, x):
+        ps = 0.0
+        if (self.saturation_pressure.type!=IncompressibleData.INCOMPRESSIBLE_NOT_SET): ps = self.psat(T, p, x)
+        if (p < 0.0): raise  ValueError("You cannot use negative pressures: {0} < {1}. ".format(p, 0.0))
+        if (p < ps) : raise ValueError("Equations are valid for liquid phase only: {0} < {1}. ".format(p, ps))
+        else        : return True
+        return False
+
+
+    def checkX(self, x):
+        if (self.xmin < 0.0 or self.xmin > 1.0): raise ValueError("Please specify the minimum concentration between 0 and 1.");
+        if (self.xmax < 0.0 or self.xmax > 1.0): raise ValueError("Please specify the maximum concentration between 0 and 1.");
+        if ((self.xmin > x) or (x > self.xmax)): raise ValueError("Your composition {0} is not between {1} and {2}.".format(x, self.xmin, self.xmax))
+        else: return True
+        return False
+    
+    def checkTPX(self, T, p, x):
+        try: 
+            return (self.checkT(T,p,x) and self.checkP(T,p,x) and self.checkX(x))
+        except ValueError as ve:
+            #print("Check failed: {0}".format(ve))
+            pass
+        return False
         
     
     def rho (self, T, p=0.0, x=0.0, c=None):
+        if not self.checkTPX(T, p, x): return np.NAN
         if c==None: 
             c=self.density.coeffs
         if self.density.type==self.density.INCOMPRESSIBLE_POLYNOMIAL:
@@ -118,6 +154,7 @@ class SolutionData(object):
         else:  raise ValueError("Unknown function.")
     
     def c   (self, T, p=0.0, x=0.0, c=None):
+        if not self.checkTPX(T, p, x): return np.NAN
         if c==None: 
             c = self.specific_heat.coeffs
         if self.specific_heat.type==self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL:
@@ -131,6 +168,7 @@ class SolutionData(object):
         return self.c(T,p,x,c)
     
     def u   (self, T, p=0.0, x=0.0, c=None):
+        if not self.checkTPX(T, p, x): return np.NAN
         if c==None: 
             c = self.specific_heat.coeffs
         if self.specific_heat.type==self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL:
@@ -143,15 +181,18 @@ class SolutionData(object):
         return self.h_u(T,p,x)
     
     def visc(self, T, p=0.0, x=0.0, c=None):
+        if not self.checkTPX(T, p, x): return np.NAN
         return self.viscosity.baseFunction(T, x, self.Tbase, self.xbase, c=c)
     
     def cond(self, T, p=0.0, x=0.0, c=None):
+        if not self.checkTPX(T, p, x): return np.NAN
         return self.conductivity.baseFunction(T, x, self.Tbase, self.xbase, c=c)
         
     def    psat(self, T, p=0.0, x=0.0, c=None):
+        if (T<=self.TminPsat): return 0.0
         return self.saturation_pressure.baseFunction(T, x, self.Tbase, self.xbase, c=c)
         
-    def Tfreeze(self, T, p=0.0, x=0.0, c=None):
+    def Tfreeze(self, T, p=0.0, x=0.0, c=None):        
         if c==None:
             c = self.T_freeze.coeffs
 
@@ -296,27 +337,31 @@ class DigitalData(SolutionData):
         if x_in!=None: # Might need update     
             if x!=None: # Both given, check if different
                 mask = np.isfinite(x)
-                if np.allclose(x[mask], x_in[mask]): 
+                if IncompressibleFitter.allClose(x[mask], x_in[mask]): 
                     if DEBUG: print("Both x-arrays are the same, no action required.")
                     updateFile = (updateFile or False) # Do not change a True value to False
                 else: 
                     updateFile = True
                     if DEBUG: print("x-arrays do not match. {0} contains \n {1} \n and will be updated with \n {2}".format(self.getFile(dataID),x,x_in))
             else: updateFile = True
-        elif x==None: raise ValueError("Could not load x from file and no x_in provided, aborting.")
+        elif x==None: 
+            if DEBUG: print("Could not load x from file {0} and no x_in provided, aborting.".format(self.getFile(dataID)))
+            return None,None,None
         else: updateFile = (updateFile or False) # Do not change a True value to False
         
         if y_in!=None: # Might need update     
             if y!=None: # Both given, check if different
                 mask = np.isfinite(y)
-                if np.allclose(y[mask], y_in[mask]): 
+                if IncompressibleFitter.allClose(y[mask], y_in[mask]): 
                     if DEBUG: print("Both y-arrays are the same, no action required.")
                     updateFile = (updateFile or False) # Do not change a True value to False
                 else: 
                     updateFile = True
                     if DEBUG: print("y-arrays do not match. {0} contains \n {1} \n and will be updated with \n {2}".format(self.getFile(dataID),y,y_in))
             else: updateFile = True
-        elif y==None: raise ValueError("Could not load y from file and no y_in provided, aborting.")
+        elif y==None:
+            if DEBUG: print("Could not load y from file {0} and no y_in provided, aborting.".format(self.getFile(dataID)))
+            return None,None,None
         else: updateFile = (updateFile or False) # Do not change a True value to False
             
         if DEBUG: print("Updating data file {0}".format(updateFile))

dev/incompressible_liquids/CPIncomp/DigitalFluids.py

@@ -10,7 +10,8 @@ parameter form.
 
 """
 from __future__ import division, print_function
-from CPIncomp.DataObjects import DigitalData  
+import numpy as np
+from CPIncomp.DataObjects import DigitalData, PureData
 
 class LiBrData(DigitalData):
     """ 
@@ -63,3 +64,198 @@ class LiBrData(DigitalData):
         self.saturation_pressure.xData,self.saturation_pressure.yData,self.saturation_pressure.data   = self.getArray(dataID=key, func=funcP, x_in=self.temperature.data, y_in=self.concentration.data,DEBUG=self.saturation_pressure.DEBUG)
         self.saturation_pressure.source = self.saturation_pressure.SOURCE_EQUATION
 
+
+
+class HyCool20(PureData,DigitalData):
+    def __init__(self):
+        DigitalData.__init__(self)
+        PureData.__init__(self)
+        
+        self.name = "HY20"
+        self.description = "HYCOOL 20, Potassium formate"
+        self.reference = "Hydro Chemicals"
+
+        self.Tmax =  50 + 273.15
+        self.Tmin = -20 + 273.15
+        self.TminPsat = self.Tmax
+        self.Tbase = 0.00 + 273.15
+        self.temperature.data         = self.getTrange()
+
+        self.density.source = self.density.SOURCE_COEFFS
+        self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
+        self.density.coeffs = np.array([[1202.2],[-0.42918]])
+
+        self.specific_heat.source = self.specific_heat.SOURCE_COEFFS
+        self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
+        self.specific_heat.coeffs = np.array([[2.955],[0.0023]])*1e3
+        
+        key = 'Cond'
+        def funcCond(T,x):
+            T = (T-self.Tbase)
+            if T <= 20: return 0.001978*T+0.5172
+            else:       return 0.001005*T+0.5368
+        self.conductivity.xData,self.conductivity.yData,self.conductivity.data = self.getArray(dataID=key,func=funcCond,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.conductivity.DEBUG)
+        self.conductivity.source = self.conductivity.SOURCE_EQUATION
+        funcCond = None
+        
+        key = 'Mu'
+        def funcMu(T,x):
+            T = (T-self.Tbase)
+            if T <= 20: return 0.07190*np.exp(524.75/(T+142.05))
+            else:       return T*(0.0005524*T - 0.06281)+2.8536
+        self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG)
+        self.viscosity.source = self.viscosity.SOURCE_EQUATION
+        funcMu = None
+
+
+class HyCool30(PureData,DigitalData):
+    def __init__(self):
+        DigitalData.__init__(self)
+        PureData.__init__(self)
+        
+        self.name = "HY30"
+        self.description = "HYCOOL 30, Potassium formate"
+        self.reference = "Hydro Chemicals"
+
+        self.Tmax =  50 + 273.15
+        self.Tmin = -30 + 273.15
+        self.TminPsat = self.Tmax
+        self.Tbase = 0.00 + 273.15
+        self.temperature.data         = self.getTrange()
+        
+        self.density.source = self.density.SOURCE_COEFFS
+        self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
+        self.density.coeffs = np.array([[1257.5],[-0.475350]])
+
+        self.specific_heat.source = self.specific_heat.SOURCE_COEFFS
+        self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
+        self.specific_heat.coeffs = np.array([[2.783],[0.0023]])*1e3
+        
+        key = 'Cond'
+        def funcCond(T,x):
+            T = (T-self.Tbase)
+            if T <= 20: return 0.001840*T+0.4980
+            else:       return 0.001000*T+0.5140
+        self.conductivity.xData,self.conductivity.yData,self.conductivity.data = self.getArray(dataID=key,func=funcCond,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.conductivity.DEBUG)
+        self.conductivity.source = self.conductivity.SOURCE_EQUATION
+        funcCond = None
+        
+        key = 'Mu'
+        def funcMu(T,x):
+            T = (T-self.Tbase)
+            if T <= 20: return 0.11100*np.exp(408.17/(T+123.15))
+            else:       return T*(0.000295*T - 0.0441)+2.6836
+        self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG)
+        self.viscosity.source = self.viscosity.SOURCE_EQUATION
+        funcMu = None
+
+
+class HyCool40(PureData,DigitalData):
+    def __init__(self):
+        DigitalData.__init__(self)
+        PureData.__init__(self)
+        
+        self.name = "HY40"
+        self.description = "HYCOOL 40, Potassium formate"
+        self.reference = "Hydro Chemicals"
+
+        self.Tmax =  20 + 273.15
+        self.Tmin = -40 + 273.15
+        self.TminPsat = self.Tmax
+        self.Tbase = 0.00 + 273.15
+        self.temperature.data         = self.getTrange()
+
+        self.density.source = self.density.SOURCE_COEFFS
+        self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
+        self.density.coeffs = np.array([[1304.5],[-0.512290]])
+
+        self.specific_heat.source = self.specific_heat.SOURCE_COEFFS
+        self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
+        self.specific_heat.coeffs = np.array([[2.646],[0.0023]])*1e3
+        
+        self.conductivity.source = self.conductivity.SOURCE_COEFFS
+        self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL
+        self.conductivity.coeffs = np.array([[0.4826],[0.001730]])
+
+        key = 'Mu'
+        def funcMu(T,x):
+            T = (T-self.Tbase)
+            return 0.07830*np.exp(498.13/(T+130.25))
+        self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG)
+        self.viscosity.source = self.viscosity.SOURCE_EQUATION
+        funcMu = None
+
+
+class HyCool45(PureData,DigitalData):
+    def __init__(self):
+        DigitalData.__init__(self)
+        PureData.__init__(self)
+        
+        self.name = "HY45"
+        self.description = "HYCOOL 45, Potassium formate"
+        self.reference = "Hydro Chemicals"
+
+        self.Tmax =  20 + 273.15
+        self.Tmin = -45 + 273.15
+        self.TminPsat = self.Tmax
+        self.Tbase = 0.00 + 273.15
+        self.temperature.data         = self.getTrange()
+        
+        self.density.source = self.density.SOURCE_COEFFS
+        self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
+        self.density.coeffs = np.array([[1328.7],[-0.530754]])
+
+        self.specific_heat.source = self.specific_heat.SOURCE_COEFFS
+        self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
+        self.specific_heat.coeffs = np.array([[2.578],[0.0023]])*1e3
+        
+        self.conductivity.source = self.conductivity.SOURCE_COEFFS
+        self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL
+        self.conductivity.coeffs = np.array([[0.4750],[0.001674]])
+
+        key = 'Mu'
+        def funcMu(T,x):
+            T = (T-self.Tbase)
+            return 0.08990*np.exp(479.09/(T+126.55))
+        self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG)
+        self.viscosity.source = self.viscosity.SOURCE_EQUATION
+        funcMu = None
+
+
+class HyCool50(PureData,DigitalData):
+    def __init__(self):
+        DigitalData.__init__(self)
+        PureData.__init__(self)
+        
+        self.name = "HY50"
+        self.description = "HYCOOL 50, Potassium formate"
+        self.reference = "Hydro Chemicals"
+
+        self.Tmax =  20 + 273.15
+        self.Tmin = -50 + 273.15
+        self.TminPsat = self.Tmax
+        self.Tbase = 0.00 + 273.15
+        self.temperature.data         = self.getTrange()
+        
+        self.density.source = self.density.SOURCE_COEFFS
+        self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
+        self.density.coeffs = np.array([[1359.0],[-0.552300]])
+
+        self.specific_heat.source = self.specific_heat.SOURCE_COEFFS
+        self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
+        self.specific_heat.coeffs = np.array([[2.498],[0.0023]])*1e3
+        
+        self.conductivity.source = self.conductivity.SOURCE_COEFFS
+        self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL
+        self.conductivity.coeffs = np.array([[0.4660],[0.001610]])
+
+        key = 'Mu'
+        def funcMu(T,x):
+            T = (T-self.Tbase)
+            res = 0.0491*np.exp(581.12/(T+129.05))
+            if T > -10: return res + 0.2
+            else:       return res  
+        self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG)
+        self.viscosity.source = self.viscosity.SOURCE_EQUATION
+        funcMu = None
+

dev/incompressible_liquids/CPIncomp/PureFluids.py

@@ -1,6 +1,6 @@
 from __future__ import division, print_function
 import numpy as np
-from CPIncomp.DataObjects import PureData  
+from CPIncomp.DataObjects import PureData  , DigitalData
    
 class TherminolD12(PureData):
     """ 
@@ -336,3 +336,4 @@ class HC10(PureData):
         self.reference   = "Dynalene data sheet"
         self.reshapeAll()
 
+    

dev/incompressible_liquids/CPIncomp/SecCoolFluids.py

@@ -45,6 +45,41 @@ class SecCoolSolutionData(DigitalData):
         self.TminPsat = self.Tmax
         
         
+        try:
+            self.density.xData,self.density.yData,self.density.data = self.getArray(dataID="Rho")
+            self.density.source = self.density.SOURCE_DATA
+        except:
+            if self.density.DEBUG: print("Could not load {}".format(self.getFile("Rho")))
+            pass
+        
+        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 
+            self.specific_heat.source = self.specific_heat.SOURCE_DATA
+        except:
+            if self.density.DEBUG: print("Could not load {}".format(self.getFile("Cp")))
+            pass
+            
+        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 
+            self.conductivity.source = self.conductivity.SOURCE_DATA
+        except:
+            if self.density.DEBUG: print("Could not load {}".format(self.getFile("Cond")))
+            pass
+        
+        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 
+            self.viscosity.source = self.viscosity.SOURCE_DATA
+        except:
+            if self.density.DEBUG: print("Could not load {}".format(self.getFile("Mu")))
+            pass
+        
+        
 #    def interp(self, tOld, xOld, dataNew):
 #        tCur = self.temperature.data
 #        xCur = self.concentration.data        
@@ -93,9 +128,7 @@ class SecCoolSolutionData(DigitalData):
         errList    = (ValueError, AttributeError, TypeError, RuntimeError)
         
         try:
-            self.density.xData,self.density.yData,self.density.data = self.getArray(dataID="Rho")
             self.density.coeffs = np.copy(std_coeffs)
-            self.density.source = self.density.SOURCE_DATA
             self.density.type   = self.density.INCOMPRESSIBLE_POLYNOMIAL
             self.density.fitCoeffs(self.Tbase,self.xbase)
         except errList as ve:
@@ -103,11 +136,7 @@ class SecCoolSolutionData(DigitalData):
             pass
         
         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 
             self.specific_heat.coeffs = np.copy(std_coeffs)
-            self.specific_heat.source = self.specific_heat.SOURCE_DATA
             self.specific_heat.type   = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
             self.specific_heat.fitCoeffs(self.Tbase,self.xbase)          
         except errList as ve:
@@ -115,84 +144,101 @@ class SecCoolSolutionData(DigitalData):
             pass 
         
         try:
-            self.conductivity.xData,self.conductivity.yData,self.conductivity.data   = self.getArray(data='Cond')
-            while np.max(self.conductivity.data[np.isfinite(self.conductivity.data)])>10: # Expect value below 1
-                self.conductivity.data *= 1e-3 
             self.conductivity.coeffs = np.copy(std_coeffs)
-            self.conductivity.source = self.conductivity.SOURCE_DATA
             self.conductivity.type   = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL
             self.conductivity.fitCoeffs(self.Tbase,self.xbase)
         except errList as ve:
             if self.conductivity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'conductivity',ve))
             pass
         
+#        try:
+#            self.viscosity.coeffs = np.copy(std_coeffs)
+#            self.viscosity.type   = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL
+#            self.viscosity.fitCoeffs(self.Tbase,self.xbase)
+#        except errList as ve:
+#            if self.viscosity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'viscosity',ve))
+#            pass
+        
         try:
-            self.viscosity.xData,self.viscosity.yData,self.viscosity.data   = self.getArray(data='Mu')
-            while np.max(self.viscosity.data[np.isfinite(self.viscosity.data)])>100: # Expect value below 10
-                self.viscosity.data *= 1e-3 
-            self.viscosity.coeffs = np.copy(std_coeffs)
-            self.viscosity.source = self.viscosity.SOURCE_DATA
-            self.viscosity.type   = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL
-            self.viscosity.fitCoeffs(self.Tbase,self.xbase)
+            tried = False
+            if len(self.viscosity.yData)==1:# and np.isfinite(fluidObject.viscosity.data).sum()<10:
+                self.viscosity.coeffs = np.array([+7e+2, -6e+1, +1e+1])
+                self.viscosity.type   = IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL
+                self.viscosity.fitCoeffs(self.Tbase,self.xbase)
+                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.type   = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL
+                self.viscosity.fitCoeffs(self.Tbase,self.xbase)
         except errList as ve:
             if self.viscosity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'viscosity',ve))
             pass
         
+        
         # reset data for getArray and read special files
         if self.xid!=self.ifrac_pure and self.xid!=self.ifrac_undefined:
             allowNegativeData_org = self.allowNegativeData
             self.allowNegativeData = True
             
-            x,_,z = self.getArray(dataID='TFreeze')
-            self.T_freeze.yData = (x - 273.15) / 100.0
-            self.T_freeze.xData = [0.0]
-            if np.min(z)<150: z += 273.15 
-            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.fitCoeffs(self.Tbase,self.xbase)
+                x,_,z = self.getArray(dataID='TFreeze')
+                self.T_freeze.yData = (x - 273.15) / 100.0
+                self.T_freeze.xData = [0.0]
+                if np.min(z)<150: z += 273.15 
+                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.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
             except errList as ve:
-                if self.T_freeze.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"T_freeze",ve))
+                if self.T_freeze.DEBUG: print("{0}: Could not load {1} data: {2}".format(self.name,"T_freeze",ve))
                 pass
 
             # reset data for getArray again
             self.allowNegativeData = allowNegativeData_org
-            x,_,z = self.getArray(dataID='Vol2Mass')            
-            massData =    z[:,0]    /100.0
-            volData  = (x - 273.15) /100.0
-            if self.xid==self.ifrac_volume:
-                _,_,self.mass2input.data = IncompressibleFitter.shapeArray(volData,axs=1)
-                self.mass2input.xData = [0.0]
-                self.mass2input.yData = massData
-                try:
-                    self.mass2input.coeffs = np.copy(std_coeffs)
-                    self.mass2input.source = self.mass2input.SOURCE_DATA
-                    self.mass2input.type   = self.mass2input.INCOMPRESSIBLE_POLYNOMIAL
-                    self.mass2input.fitCoeffs(self.Tbase,self.xbase)
-                except errList as ve:
-                    if self.mass2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"mass2input",ve))
-                    pass
-            elif self.xid==self.ifrac_mass:
-                _,_,self.volume2input.data = IncompressibleFitter.shapeArray(massData,axs=1)
-                self.volume2input.xData = [0.0]
-                self.volume2input.yData = volData
-                try:
-                    self.volume2input.coeffs = np.copy(std_coeffs)
-                    self.volume2input.source = self.volume2input.SOURCE_DATA
-                    self.volume2input.type   = self.volume2input.INCOMPRESSIBLE_POLYNOMIAL
-                    self.volume2input.fitCoeffs(self.Tbase,self.xbase)
-                except errList as ve:
-                    if self.volume2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"volume2input",ve))
-                    pass
-            else:
-                raise ValueError("Unknown xid specified.")
+            try:
+                x,_,z = self.getArray(dataID='Vol2Mass')            
+                massData =    z[:,0]    /100.0
+                volData  = (x - 273.15) /100.0
+                
+                if self.xid==self.ifrac_volume:
+                    _,_,self.mass2input.data = IncompressibleFitter.shapeArray(volData,axs=1)
+                    self.mass2input.xData = [0.0]
+                    self.mass2input.yData = massData
+                    try:
+                        self.mass2input.coeffs = np.copy(std_coeffs)
+                        self.mass2input.source = self.mass2input.SOURCE_DATA
+                        self.mass2input.type   = self.mass2input.INCOMPRESSIBLE_POLYNOMIAL
+                        self.mass2input.fitCoeffs(self.Tbase,self.xbase)
+                    except errList as ve:
+                        if self.mass2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"mass2input",ve))
+                        pass
+                elif self.xid==self.ifrac_mass:
+                    _,_,self.volume2input.data = IncompressibleFitter.shapeArray(massData,axs=1)
+                    self.volume2input.xData = [0.0]
+                    self.volume2input.yData = volData
+                    try:
+                        self.volume2input.coeffs = np.copy(std_coeffs)
+                        self.volume2input.source = self.volume2input.SOURCE_DATA
+                        self.volume2input.type   = self.volume2input.INCOMPRESSIBLE_POLYNOMIAL
+                        self.volume2input.fitCoeffs(self.Tbase,self.xbase)
+                    except errList as ve:
+                        if self.volume2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"volume2input",ve))
+                        pass
+                else:
+                    raise ValueError("Unknown xid specified.")
+            except errList as ve:
+                if self.mass2input.DEBUG or self.volume2input.DEBUG: print("{0}: Could not load {1} data: {2}".format(self.name,"Vol2Mass",ve))
+                pass
 
         
     # Redefine some functions to avoid data loss
@@ -331,6 +377,13 @@ class SecCoolSolutionData(DigitalData):
         sec += [SecCoolSolutionData(sFile='Dowtherm Q'   ,sFolder='xPure',name='DowQ2',desc='Dowtherm Q, Diphenylethane/alkylated aromatics' ,ref='Dow Chemicals, SecCool software')]
         print(", {0}".format(sec[-1].name), end="")
         
+        sec += [SecCoolIceData(sFile='IceEA'   ,sFolder='xMass',name='IceEA',desc='Ice slurry with Ethanol' ,ref='Danish Technological Institute, SecCool software')]
+        print(", {0}".format(sec[-1].name), end="")
+        sec += [SecCoolIceData(sFile='IceNA'   ,sFolder='xMass',name='IceNA',desc='Ice slurry with NaCl' ,ref='Danish Technological Institute, SecCool software')]
+        print(", {0}".format(sec[-1].name), end="")
+        sec += [SecCoolIceData(sFile='IcePG'   ,sFolder='xMass',name='IcePG',desc='Ice slurry with Propylene Glycol' ,ref='Danish Technological Institute, SecCool software')]
+        print(", {0}".format(sec[-1].name), end="")
+        
         sec += [ThermogenVP1869()]
         print(", {0}".format(sec[-1].name), end="")
         sec += [Freezium()]
@@ -346,13 +399,107 @@ class SecCoolSolutionData(DigitalData):
         sec += [AS55()]
         print(", {0}".format(sec[-1].name), end="")
         
-        
-        
         print(" ... done")
         
         return sec
 
 
+
+class SecCoolIceData(SecCoolSolutionData):
+    """ 
+    A base class that can be fed with a fluid ID from SecCool
+    to read data files sitting in data/SecCool/xTables.
+    """ 
+    def __init__(self,sFile=None,sFolder=None,name=None,desc=None,ref='Danish Technological Institute, SecCool software'):
+        SecCoolSolutionData.__init__(self,sFile=sFile,sFolder=sFolder,name=name,desc=desc,ref=ref)
+        
+        #self.density.xData,self.density.yData,self.density.data = self.getArray(dataID="Rho")
+        #self.density.source = self.density.SOURCE_DATA
+        
+        self.specific_heat.xData,self.specific_heat.yData,self.specific_heat.data = self.getArray(dataID='Hfusion')
+        self.specific_heat.source = self.specific_heat.SOURCE_DATA
+        
+        #self.conductivity.xData,self.conductivity.yData,self.conductivity.data   = self.getArray(dataID='Cond')
+        #self.conductivity.source = self.conductivity.SOURCE_DATA
+        
+        #self.viscosity.xData,self.viscosity.yData,self.viscosity.data   = self.getArray(dataID='Mu')
+        #self.viscosity.source = self.viscosity.SOURCE_DATA
+        
+        
+#    def fitFluid(self):
+#        if self.Tbase==None:
+#            self.Tbase = (self.Tmin + self.Tmax) / 2.0
+#        if self.xbase==None:
+#            self.xbase = (self.xmin + self.xmax) / 2.0
+#        
+#        std_coeffs = np.zeros((4,6))
+#        errList    = (ValueError, AttributeError, TypeError, RuntimeError)
+#        
+#        try:
+#            self.density.coeffs = np.copy(std_coeffs)
+#            self.density.type   = self.density.INCOMPRESSIBLE_POLYNOMIAL
+#            self.density.fitCoeffs(self.Tbase,self.xbase)
+#        except errList as ve:
+#            if self.density.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'density',ve))
+#            pass
+#        
+#        try:
+#            self.specific_heat.coeffs = np.copy(std_coeffs)
+#            self.specific_heat.type   = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
+#            self.specific_heat.fitCoeffs(self.Tbase,self.xbase)          
+#        except errList as ve:
+#            if self.specific_heat.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'specific heat',ve))
+#            pass 
+#        
+#        try:
+#            self.conductivity.coeffs = np.copy(std_coeffs)
+#            self.conductivity.type   = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL
+#            self.conductivity.fitCoeffs(self.Tbase,self.xbase)
+#        except errList as ve:
+#            if self.conductivity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'conductivity',ve))
+#            pass
+#        
+#        try:
+#            self.viscosity.coeffs = np.copy(std_coeffs)
+#            self.viscosity.type   = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL
+#            self.viscosity.fitCoeffs(self.Tbase,self.xbase)
+#        except errList as ve:
+#            if self.viscosity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'viscosity',ve))
+#            pass
+
+        
+    # Redefine some functions to avoid data loss
+    def getFile(self, data):
+        return os.path.join(os.path.dirname(__file__), 'data','SecCool','xTables', self.sFolder, self.sFile+"_"+data+".csv")
+        
+    
+    def getFromFile(self, data):
+        fullPath = self.getFile(data)
+        content = np.loadtxt(fullPath,dtype=type('string'),delimiter=",",skiprows=3)
+        r,c,res = IncompressibleFitter.shapeArray(content)
+#        numbers = res.astype(np.float)
+        numbers = np.zeros((r,c))
+        for i in range(r):
+            for j in range(c):
+                nu = np.NAN
+                try:
+                    nu = np.float(res[i,j])
+                    if i==0: nu *= 1e-2 # Percent to fraction
+                    if not self.allowNegativeData and nu<0:
+                        nu = np.NAN # invalid entries
+                except (ValueError, TypeError) as ve:
+                    if False: print("Could not convert entry: {0}".format(ve))
+                    if i==0: nu = 0.0 # Dummy for tables without concentration (TFreeze and Vol2Mass)
+                    pass
+                numbers[i,j] = nu
+        if numbers[1,0]>numbers[-1,0]: 
+            numbers[1:,:] = numbers[1:,:][::-1,:]
+        if numbers[0,1]>numbers[0,-1]: 
+            numbers[:,1:] = numbers[:,1:][:,::-1]
+        return numbers 
+
+
+
 class ThermogenVP1869(PureData,DigitalData):
     """ 
     Source: SecCool Software
@@ -731,361 +878,7 @@ class AS55(PureData,DigitalData):
 
 
 
-    
-    
-if __name__ == '__main__':   
-    
-    print("The main method converts the AspenTemper coefficients from SecCool")
-    names = []
-    descs = []
-    refs  = []
-    
-    Tmax = []
-    Tmin = []
-    TminPsat = []
-    Tbase = []
-    
-    cond = []
-    cp   = []
-    mu   = []
-    rho  = []
-    
-    ID = 10
-    names.append("AS{0}".format(ID))
-    descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID))
-    refs.append ("Aspen Petroleum AB, SecCool software")
-    
-    Tmax.append("30 + 273.15")
-    Tmin.append("-{0}".format(ID))
-    TminPsat.append("self.Tmax")
-    Tbase.append("0.00 + 273.15")
-    
-    cond.append ("0.001483*(273+T)+0.1094")
-    cp.append   ("1000*(3.54183+T*(0.00201-0.0000132589*T))")
-    mu.append   ("2.80154921+T*(-0.10631376+T*(0.00235381-0.0000211481*T))")
-    rho.append  ("1090+T*(-0.2+T*(1.66533E-16-1.89735E-18*T))")
-    
-    ID = 20
-    names.append("AS{0}".format(ID))
-    descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID))
-    refs.append ("Aspen Petroleum AB, SecCool software")
-        
-    Tmax.append("30 + 273.15")
-    Tmin.append("-{0}".format(ID))
-    TminPsat.append("self.Tmax")
-    Tbase.append("0.00 + 273.15")
-    
-    cond.append ("0.001342*(273+T)+0.1144")
-    cp.append   ("1000*(3.26252+T*(0.00286-0.0000122673*T))")
-    mu.append   ("(0.97244+(7.14199*exp((-20-T)/18.6014)))")
-    rho.append  ("1147+T*(-0.22142857+T*(-0.00142857+2.98156E-19*T))")
 
-
-    ID = 30
-    names.append("AS{0}".format(ID))
-    descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID))
-    refs.append ("Aspen Petroleum AB, SecCool software")
-        
-    Tmax.append("30 + 273.15")
-    Tmin.append("-{0}".format(ID))
-    TminPsat.append("self.Tmax")
-    Tbase.append("0.00 + 273.15")
-    
-    cond.append ("0.001256*(273+T)+0.1175")
-    cp.append   ("1000*(3.07504+T*(0.00299-0.0000270232*T))")
-    mu.append   ("(1.30143+(16.01368*exp((-30-T)/16.69989)))")
-    rho.append  ("1183.85930736+T*(-0.31085859+T*(-0.00103896+0.0000176768*T))")
-
-
-    ID = 40
-    names.append("AS{0}".format(ID))
-    descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID))
-    refs.append ("Aspen Petroleum AB, SecCool software")
-        
-    Tmax.append("30 + 273.15")
-    Tmin.append("-{0}".format(ID))
-    
-    cond.append ("0.001099*(273+T)+0.1433")
-    cp.append   ("1000*(2.97788+T*(0.00228-0.0000387227*T))")
-    mu.append   ("(39.11536*exp((-40-T)/9.99495)+(12.41564*exp((-40-T)/38.45577)))")
-    rho.append  ("1214.83982684+T*(-0.37819865+T*(-0.00109307+0.000016835*T))")
-
-    ID = 55
-    names.append("AS{0}".format(ID))
-    descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID))
-    refs.append ("Aspen Petroleum AB, SecCool software")
-        
-    Tmax.append("30 + 273.15")
-    Tmin.append("-{0}".format(ID))
-    TminPsat.append("self.Tmax")
-    Tbase.append("0.00 + 273.15")
-    
-    cond.append ("(273+T)*(0.000002287*(273+T)-0.0003108)+0.3402")
-    cp.append   ("1000*(2.83985+T*(0.00229-0.0000248618*T))")
-    mu.append   ("(317.40673*exp((-55-T)/7.24125)+(51.22151*exp((-55-T)/26.28052)))")
-    rho.append  ("1249.7534665+T*(-0.47629615+T*(-0.00117189+0.0000198824*T))")
-    
-    spa = "        "
-    
-    from sympy import symbols, expand
-    x, T = symbols('x T') 
-    
-    for i in range(len(names)):
-        print("class {0}(PureData):\n    def __init__(self):".format(names[i]))
-        print("{0}PureData.__init__(self)".format(spa))
-        print("{0}self.name = \"{1}\"".format(spa,names[i]))
-        print("{0}self.description = \"{1}\"".format(spa,descs[i]))
-        print("{0}self.reference = \"{1}\"".format(spa,refs[i]))
-        print("")
-        print("{0}self.Tmax = {1}".format(spa,Tmax[i]))
-        print("{0}self.Tmin = {1}".format(spa,Tmin[i]))
-        print("{0}self.TminPsat = self.Tmax".format(spa))
-        print("{0}self.Tbase = 0.00 + 273.15".format(spa))
-        print("")
-        print("{0}self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL".format(spa))
-        print("{0}self.density.coeffs = np.array([{1}])[::-1]".format(spa,expand(rho[i])))
-        print("")
-        print("{0}self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL".format(spa))
-        print("{0}self.specific_heat.coeffs = np.array([{1}])[::-1]".format(spa,expand(cp[i])))
-        print("")
-        print("{0}self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL".format(spa))
-        print("{0}self.conductivity.coeffs = np.array([{1}])[::-1]".format(spa,expand(cond[i])))
-        print("")
-        print("{0}self.viscosity.type = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL".format(spa))
-        print("{0}self.viscosity.coeffs = np.array([{1}])[::-1]".format(spa,expand(mu[i])))
-        print("")
-        print("{0}def fitFluid(self):".format("    "))
-        print("{0}key = \"Mu\"".format(spa))
-        print("{0}def funcMu(T,x):".format(spa))
-        print("")
-        print("")
-        print("{0}funcMu=None".format(spa))
-        print("")
-        print("")
-        print("")
-        print("")
-        
-#        PureData.__init__(self)
-#        self.name = "AS10"
-#        self.description = "Aspen Temper -10, Potassium acetate/formate"
-#        self.reference = "Aspen Petroleum AB, SecCool software"
-#        
-#        self.Tmax =  30 + 273.15
-#        self.Tmin = -10 + 273.15
-#        self.TminPsat =  self.Tmax
-#    
-#        self.Tbase =   0.00 + 273.15
-#       
-#        self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
-#        self.density.coeffs = np.array([[945.5454545],[-1.054545455]])
-#        
-#        self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL
-#        self.specific_heat.coeffs = np.array([[2.322218182],[0.003843636]])*1000
-#
-#        self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL
-#        self.conductivity.coeffs = np.array([[0.001483*273.15+0.1094],[0.001483]])
-#        
-#        self.temperature.data         = self.getTrange()
-#        self.concentration.data       = np.array([0]) # mass fraction
-
-
-#
-#
-#{ THyCool_20 }
-#
-#constructor THyCool_20.Create;
-#begin
-#  inherited;
-#  TFreeze   := -20;
-#  FluidType := 'Potassium Formate';
-#  TradeName := 'HYCOOL 20';
-#  Reference := 'Hydro Chemicals';
-#  TMin      := -20;
-#  TMax      := 50;
-#end;
-#
-#function THyCool_20.GetCondT(T: Double): Double;
-#begin
-#  if T <= 20 then
-#    Result := 0.001978*T+0.5172
-#  else
-#    Result := 0.001005*T+0.5368;
-#end;
-#
-#function THyCool_20.GetCpT(T: Double): Double;
-#begin
-#  Result := 1000*(0.0023*T+2.955);
-#end;
-#
-#function THyCool_20.GetMuT(T: Double): Double;
-#begin
-#  if T <= 20 then
-#    Result := 0.07190*exp(524.75/(T+142.05))
-#  else
-#    Result := T*(0.0005524*T - 0.06281)+2.8536;
-#end;
-#
-#function THyCool_20.GetRhoT(T: Double): Double;
-#begin
-#  Result := -0.429180*T+1202.2;
-#end;
-#
-#{ THyCool_30 }
-#
-#constructor THyCool_30.Create;
-#begin
-#  inherited;
-#  TFreeze   := -30;
-#  FluidType := 'Potassium Formate';
-#  TradeName := 'HYCOOL 30';
-#  Reference := 'Hydro Chemicals';
-#  TMin      := -30;
-#  TMax      := 50;
-#end;
-#
-#function THyCool_30.GetCondT(T: Double): Double;
-#begin
-#  if T <= 20 then
-#    Result := 0.001840*T+0.4980
-#  else
-#    Result := 0.001000*T+0.514;
-#end;
-#
-#function THyCool_30.GetCpT(T: Double): Double;
-#begin
-#  Result := 1000*(0.0023*T+2.783);
-#end;
-#
-#function THyCool_30.GetMuT(T: Double): Double;
-#begin
-#  if T <= 20 then
-#    Result := 0.11100*exp(408.17/(T+123.15))
-#  else
-#    Result := T*(0.000295*T - 0.0441)+2.6836;
-#end;
-#
-#function THyCool_30.GetRhoT(T: Double): Double;
-#begin
-#  Result := -0.475350*T+1257.5;
-#end;
-#
-#{ THyCool_40 }
-#
-#constructor THyCool_40.Create;
-#begin
-#  inherited;
-#  TFreeze   := -40;
-#  FluidType := 'Potassium Formate';
-#  TradeName := 'HYCOOL 40';
-#  Reference := 'Hydro Chemicals';
-#  TMin      := -40;
-#  TMax      := 20;
-#end;
-#
-#function THyCool_40.GetCondT(T: Double): Double;
-#begin
-#  Result := 0.001730*T+0.4826;
-#end;
-#
-#function THyCool_40.GetCpT(T: Double): Double;
-#begin
-#  Result := 1000*(0.0023*T+2.646);
-#end;
-#
-#function THyCool_40.GetMuT(T: Double): Double;
-#begin
-#  Result := 0.07830*exp(498.13/(T+130.25));
-#end;
-#
-#function THyCool_40.GetRhoT(T: Double): Double;
-#begin
-#  Result := -0.512290*T+1304.5;
-#end;
-#
-#{ THyCool_45 }
-#
-#constructor THyCool_45.Create;
-#begin
-#  inherited;
-#  TFreeze   := -45;
-#  FluidType := 'Potassium Formate';
-#  TradeName := 'HYCOOL 45';
-#  Reference := 'Hydro Chemicals';
-#  TMin      := -40;
-#  TMax      := 20;
-#end;
-#
-#function THyCool_45.GetCondT(T: Double): Double;
-#begin
-#  Result := 0.001674*T+0.4750;
-#end;
-#
-#function THyCool_45.GetCpT(T: Double): Double;
-#begin
-#  Result := 1000*(0.0023*T+2.578);
-#end;
-#
-#function THyCool_45.GetMuT(T: Double): Double;
-#begin
-#  Result := 0.08990*exp(479.09/(T+126.55));
-#end;
-#
-#function THyCool_45.GetRhoT(T: Double): Double;
-#begin
-#  Result := -0.530754*T+1328.7;
-#end;
-#
-#{ THyCool_50 }
-#
-#constructor THyCool_50.Create;
-#begin
-#  inherited;
-#  TFreeze   := -50;
-#  FluidType := 'Potassium Formate';
-#  TradeName := 'HYCOOL 50';
-#  Reference := 'Hydro Chemicals';
-#  TMin      := -50;
-#  TMax      := 20;
-#end;
-#
-#function THyCool_50.GetCondT(T: Double): Double;
-#begin
-#  Result := 0.001610*T+0.4660;
-#end;
-#
-#function THyCool_50.GetCpT(T: Double): Double;
-#begin
-#  Result := 1000*(0.0023*T+2.498);
-#end;
-#
-#function THyCool_50.GetMuT(T: Double): Double;
-#begin
-#  Result := 0.0491*exp(581.12/(T+129.05));
-#  if T > -10 then
-#    Result := Result+0.2;
-#end;
-#
-#function THyCool_50.GetRhoT(T: Double): Double;
-#begin
-#  Result := -0.552300*T+1359.0;
-#end;
-#
-#        
-#        
-#        
-#        
-#        
-#        
-#        
-#        
-#        
-#        
-#        
-#        
-#        
-        
-        
-        
         
         
         

dev/incompressible_liquids/CPIncomp/WriterObjects.py

@@ -5,8 +5,9 @@ import hashlib, os, json, sys
 import matplotlib.pyplot as plt 
 import matplotlib.gridspec as gridspec
 from CPIncomp.DataObjects import SolutionData
-from CPIncomp.BaseObjects import IncompressibleData
+from CPIncomp.BaseObjects import IncompressibleData, IncompressibleFitter
 from matplotlib.patches import Rectangle
+from matplotlib.ticker import MaxNLocator
 
 class SolutionDataWriter(object):
     """ 
@@ -20,8 +21,6 @@ class SolutionDataWriter(object):
     
     def fitAll(self, fluidObject=SolutionData()):
         
-        
-        
         if fluidObject.Tbase==None:
             fluidObject.Tbase = (fluidObject.Tmin + fluidObject.Tmax) / 2.0
             
@@ -71,10 +70,10 @@ class SolutionDataWriter(object):
             fluidObject.viscosity.setxyData(tData,xData)
             tried = False
             if len(fluidObject.viscosity.yData)==1:# and np.isfinite(fluidObject.viscosity.data).sum()<10:
-                fluidObject.viscosity.coeffs = np.array([+7e+2, -6e+1, +1e+1])
+                fluidObject.viscosity.coeffs = np.array([+5e+2, -6e+1, +1e+1])
                 fluidObject.viscosity.type   = IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL
                 fluidObject.viscosity.fitCoeffs(tBase,xBase)
-                if fluidObject.viscosity.coeffs==None or np.allclose(fluidObject.viscosity.coeffs, np.array([+7e+2, -6e+1, +1e+1])): # Fit failed
+                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))
@@ -86,13 +85,17 @@ class SolutionDataWriter(object):
         
         try:
             fluidObject.saturation_pressure.setxyData(tData,xData)
+            tried = False
             if len(fluidObject.saturation_pressure.yData)==1:# and np.isfinite(fluidObject.saturation_pressure.data).sum()<10:
                 fluidObject.saturation_pressure.coeffs = np.array([-5e+3, +6e+1, -1e+1]) 
                 fluidObject.saturation_pressure.type   = IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL
-            else:
+                fluidObject.saturation_pressure.fitCoeffs(tBase,xBase)
+                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.type   = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL
-            fluidObject.saturation_pressure.fitCoeffs(tBase,xBase)
+                fluidObject.saturation_pressure.fitCoeffs(tBase,xBase)
         except errList as ve:
             if fluidObject.saturation_pressure.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(fluidObject.name,'saturation pressure',ve))
             pass
@@ -307,7 +310,7 @@ class SolutionDataWriter(object):
         
         pos = np.isfinite(B_a)
         pos2 = (B_a>eps)
-        result = np.zeros_like(A_a)
+        result = np.ones_like(A_a)*np.NAN
         
         result[pos & pos2] = (A_a[pos & pos2]-B_a[pos & pos2])/B_a[pos & pos2]
 
@@ -416,9 +419,7 @@ class SolutionDataWriter(object):
                         zError[i,j]= func(tData[i],pData,xData[j])
                                                    
                 zError = self.relError(zData, zError) * 1e2
-                
-                if xFunction: axis = 1
-                else: axis = 0
+
                   
                 ## Find the column with the largest single error
                 #maxVal = np.amax(zError, axis=0) # largest error per column
@@ -426,17 +427,21 @@ class SolutionDataWriter(object):
                 ## Find the column with the largest total error
                 #totVal = np.sum(zError, axis=0) # summed error per column
                 #col2plot = np.argmax(totVal) # largest error in row
-                # Find the column with the largest average error
-                avgVal = np.average(zError, axis=axis) # summed error per column
-                set2plot = np.argmax(avgVal) # largest error in row
+                # Find the column with the largest average error                
                 if xFunction: 
+                    #avgVal = np.average(zError, axis=1) # summed error per column
+                    #set2plot = np.argmax(avgVal) # largest error in row
+                    set2plot = int(np.round(r/2.0))
                     tData = np.array([tData[set2plot]])
                     zData = zData[set2plot]
                     zError= zError[set2plot]
                 else:
+                    #avgVal = np.average(zError, axis=0) # summed error per column
+                    #set2plot = np.argmax(avgVal) # largest error in row
+                    set2plot = int(np.round(c/2.0))
                     xData = np.array([xData[set2plot]])
                     zData = zData.T[set2plot]
-                    zError= zError.T[set2plot]                
+                    zError= zError.T[set2plot]
             else:
                 raise ValueError("You have to provide data and a fitted function.")
         
@@ -472,6 +477,10 @@ class SolutionDataWriter(object):
         xFunc = xData
         pFunc = pData
         zFunc = None
+        zMiMa = None
+        xFree = xData
+        tFree = tData
+        zFree = None
         
         if func!=None:
             if len(tFunc)<points and len(tFunc)>1:
@@ -488,7 +497,30 @@ class SolutionDataWriter(object):
                 zFunc = np.array(zFunc.flat)
             else:
                 raise ValueError("Cannot plot non-flat arrays!")
-        
+            
+            if xFunction:
+                tMiMa = np.array([solObj.Tmin, solObj.Tmax])
+                xMiMa = xFunc
+            else:
+                tMiMa = tFunc
+                xMiMa = np.array([solObj.xmin, solObj.xmax])
+                
+            zMiMa = np.zeros((len(tMiMa),len(xMiMa)))
+            for i in range(len(tMiMa)):
+                for j in range(len(xMiMa)):
+                    zMiMa[i,j] = func(tMiMa[i],pFunc,xMiMa[j])
+                    
+            if not xFunction: # add the freezing front
+                if solObj.T_freeze.type!=IncompressibleData.INCOMPRESSIBLE_NOT_SET:
+                    cols = len(tMiMa)
+                    conc = np.linspace(solObj.xmin, solObj.xmax, num=cols)
+                    tFree = np.zeros_like(conc)
+                    zFree = np.zeros_like(conc)
+                    for i in range(cols):
+                        tFree[i] = solObj.Tfreeze(10.0, p=pFunc, x=conc[i])
+                        zFree[i] = func(tFree[i],pFunc,conc[i])
+                    #zMiMa = np.hstack((zMiMa,temp.reshape((len(conc),1)))) 
+
         fitFormatter = {}
         fitFormatter['color'] = 'red'
         fitFormatter['ls'] = 'solid'
@@ -500,12 +532,20 @@ class SolutionDataWriter(object):
         errorFormatter['alpha'] = 0.25
         
         pData = None
+        pFree = None
         if xFunction:
             pData = xData
             pFunc = xFunc
+            pMiMa = xMiMa
+            zMiMa = zMiMa.T
+            #pFree = xFree
+            #zFree = zFree.T
         else:
             pData = tData - 273.15 
-            pFunc = tFunc - 273.15 
+            pFunc = tFunc - 273.15
+            pMiMa = tMiMa - 273.15  
+            #zMiMa = zMiMa
+            pFree = tFree - 273.15
             
         if zData!=None and axVal!=None: 
             axVal.plot(pData, zData, label='data', **dataFormatter)
@@ -513,6 +553,12 @@ class SolutionDataWriter(object):
         if zFunc!=None and axVal!=None:
             axVal.plot(pFunc, zFunc, label='function' , **fitFormatter)
             
+        if zMiMa!=None and axVal!=None:
+            axVal.plot(pMiMa, zMiMa, alpha=0.25, ls=':', color=fitFormatter["color"])
+            
+        if zFree!=None and axVal!=None:
+            axVal.plot(pFree, zFree, alpha=0.25, ls=':', color=fitFormatter["color"])
+            
         if zError!=None and axErr!=None:
             axErr.plot(pData, zError, label='error' , **errorFormatter)
             if solObj.xid!=solObj.ifrac_pure and not xFunction:
@@ -671,7 +717,7 @@ class SolutionDataWriter(object):
         density_axis.set_ylabel(r'Density [$\mathdefault{kg/m^3\!}$]')
         density_axis.set_xlabel(tempLabel)
         density_error.set_ylabel(errLabel)
-        if solObj.density.source!=solObj.density.SOURCE_NOT_SET:        
+        if solObj.density.source!=solObj.density.SOURCE_NOT_SET:
             self.plotValues(density_axis,density_error,solObj=solObj,dataObj=solObj.density,func=solObj.rho)
         else:
             raise ValueError("Density data has to be provided!")
@@ -751,13 +797,15 @@ class SolutionDataWriter(object):
         #mass2input_axis   = plt.subplot2grid((3,2), (2,0))
         #volume2input_axis = plt.subplot2grid((3,2), (2,0))
         
-        # Set a minimum error level
+        # Set a minimum error level and do some more formatting
         minAbsErrorScale = 0.05 # in per cent
         for a in fig.axes:
             if a.get_ylabel()==errLabel:
                 mi,ma = a.get_ylim()
                 if mi>-minAbsErrorScale: a.set_ylim(bottom=-minAbsErrorScale)
                 if ma< minAbsErrorScale: a.set_ylim(   top= minAbsErrorScale)
+            a.xaxis.set_major_locator(MaxNLocator(5))
+            #a.yaxis.set_major_locator(MaxNLocator(7))
                 
         
         
@@ -779,7 +827,7 @@ class SolutionDataWriter(object):
         
         table_axis.legend( 
           legVal, legKey, 
-          bbox_to_anchor=(0.0, 0.015, 1., 0.015), 
+          bbox_to_anchor=(0.0, -0.025, 1., -0.025), 
           ncol=len(legKey), mode="expand", borderaxespad=0.,
           numpoints=1)
         #table_axis.legend(handles, labels, bbox_to_anchor=(0.0, -0.1), loc=2, ncol=3)

dev/incompressible_liquids/all_incompressibles.py

@@ -10,22 +10,28 @@ if __name__ == '__main__':
     writer = SolutionDataWriter()
     doneObjs = []
     
-#    # To debug single fluids
-#    from CPIncomp.SecCoolFluids import Freezium
-#    #solObjs = [SecCoolSolutionData(sFile='Melinder, Ammonia'            ,sFolder='xMass',name='MAM2',desc='Melinder, Ammonia'            ,ref='Melinder-BOOK-2010, 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)
-#    #
-#    ##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
+    #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)
     
     fluidObjs = getExampleNames(obj=True)
     examplesToFit = ["ExamplePure","ExampleSolution","ExampleDigital","ExampleDigitalPure"]

src/Backends/Incompressible/IncompressibleFluid.cpp

@@ -439,26 +439,21 @@ double IncompressibleFluid::T_u   (double Umass, double p, double x){
  *  freezing point calculation. This is not necessarily
  *  defined for all fluids, default values do not cause errors. */
 bool IncompressibleFluid::checkT(double T, double p, double x) {
-	if (Tmin <= 0.) {
-		throw ValueError("Please specify the minimum temperature.");
-	} else if (Tmax <= 0.) {
-		throw ValueError("Please specify the maximum temperature.");
-	} else if ((Tmin > T) || (T > Tmax)) {
-		throw ValueError(
+	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));
+	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 {
-		double TF = 0.0;
-		if (T_freeze.type!=IncompressibleData::INCOMPRESSIBLE_NOT_SET) {
-			TF = Tfreeze(p, x);
-		}
-		if ( T<Tfreeze(p, x)) {
-			throw ValueError(
-				format("Your temperature %f is below the freezing point of %f.",
-						T, Tfreeze(p, x)));
-		} else {
-			return true;
-		}
+		return true;
 	}
 	return false;
 }