In [1]:
import os, sys, inspect

currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
if os.path.basename(currentdir) == 'notebooks':
    parentdir = os.path.dirname(currentdir)
    sys.path.insert(0,parentdir)
elif os.path.basename(currentdir) == 'celldeposit-condenser':
    currentdir = os.path.join(currentdir, 'notebooks')
    os.chdir(currentdir)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from water_at_saturation_properties import density, enthalpy, heat_capacity, conductivity, viscosity
from water_at_saturation_properties import saturation_temperature, vapour_density, vaporization_enthalpy

In [2]:
def simplified_model_for_fit(x, a00, a01, a02, a03):
    
    c = np.array([a00, a01, a02, a03,])
    
    return simplified_model(x, c)

In [3]:
def simplified_model(x,c):

    f = 0.

    for i in range(c.shape[0]):
        f += c[i] * (x ** i)

    return f

In [4]:
def calculate_output(func,X,*args):
    Z = np.zeros_like(X)
    for i in range(X.shape[0]):                     
        Z[i,] = func(X[i],*args)
    return Z

In [5]:
def fit(model, X):

    Z = calculate_output(model,X)
    c, pcov = curve_fit(simplified_model_for_fit, X, Z)
    
    Z_fit = simplified_model(X,c)
    
    # residual sum of squares
    ss_res = np.sum((Z - Z_fit) ** 2)

    # total sum of squares
    ss_tot = np.sum((Z - np.mean(Z)) ** 2)

    # r-squared
    r2 = 1 - (ss_res / ss_tot)
    
    return (c, Z_fit, Z, r2)

In [6]:
from iapws.iapws97 import _Region4, _Region1
from iapws._iapws import _ThCond, _Viscosity

In [7]:
limits = [-1, 0]
npoints = 25
X = 1e5*np.logspace(limits[0], limits[1], num=npoints)
print(X)
out = fit(vapour_density, X)
print(out)


[ 10000.          11006.94171252  12115.27658629  13335.21432163
  14677.99267622  16155.9809844   17782.79410039  19573.41781488
  21544.34690032  23713.73705662  26101.57215683  28729.84833354
  31622.77660168  34807.00588428  38311.86849557  42169.65034286
  46415.88833613  51089.69774507  56234.13251903  61896.58188913
  68129.2069058   74989.42093325  82540.4185268   90851.75756517
 100000.        ]
(array([ 4.62513512e-03,  6.47401547e-06, -9.47682219e-12,  3.32572527e-17]), array([0.06845086, 0.07478045, 0.08172776, 0.08935114, 0.09771414,
       0.10688585, 0.1169414 , 0.12796239, 0.14003739, 0.15326252,
       0.16774204, 0.18358906, 0.20092634, 0.21988724, 0.24061688,
       0.26327361, 0.28803084, 0.31507956, 0.3446315 , 0.37692362,
       0.41222403, 0.45084019, 0.49313017, 0.53951815, 0.59051571]), array([0.06816373, 0.07460199, 0.08164584, 0.08935202, 0.09778257,
       0.10700532, 0.11709447, 0.12813112, 0.14020397, 0.15341   ,
       0.16785529, 0.18365582, 0.20093843, 0.21984186, 0.24051783,
       0.2631323 , 0.28786683, 0.31492003, 0.34450918, 0.37687205,
       0.41226884, 0.45098433, 0.49333024, 0.53964787, 0.59031092]), 0.9999991773107247)

In [8]:
print("VAPOUR DENSITY")
c, Z_fit, Z, r2 = fit(vapour_density, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict = dict()
c_dict['vapour_density'] = c


VAPOUR DENSITY
Parameters: [ 4.62513512e-03  6.47401547e-06 -9.47682219e-12  3.32572527e-17]
R2: 0.9999991773107247

In [9]:
print("VAPOUR TOTAL COMPRESSIBILITY")
c_v_density = c_dict['vapour_density']
c = np.zeros((c_v_density.shape[0]-1,))
i = 0
for ci in c_v_density:
    if i > 0:
        c[i-1] = i*ci
    i += 1
Z_fit = simplified_model(X,c)
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c_v_density)
print("Parameters:",c)
print("R2:", r2)
c_dict['vapour_total_compressibility'] = c


VAPOUR TOTAL COMPRESSIBILITY
Parameters: [ 4.62513512e-03  6.47401547e-06 -9.47682219e-12  3.32572527e-17]
Parameters: [ 6.47401547e-06 -1.89536444e-11  9.97717582e-17]
R2: 0.9999991773107247

In [10]:
print("SATURATION TEMPERATURE")
c, Z_fit, Z, r2 = fit(saturation_temperature, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict['saturation_temperature'] = c


SATURATION TEMPERATURE
Parameters: [ 3.04528132e+02  1.74387838e-03 -1.88042391e-08  8.25850104e-14]
R2: 0.9989376467549783

In [11]:
print("DENSITY")
c, Z_fit, Z, r2 = fit(density, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict['density'] = c


DENSITY
Parameters: [ 9.96906358e+02 -8.11519963e-04  7.28615211e-09 -3.02080951e-14]
R2: 0.9996968680211483

In [12]:
print("ENTHALPY")
c, Z_fit, Z, r2 = fit(enthalpy, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict['enthalpy'] = c


ENTHALPY
Parameters: [ 1.31488947e+05  7.28807358e+00 -7.84449119e-05  3.44521411e-10]
R2: 0.9989472814741913

In [13]:
print("VISCOSITY")
c, Z_fit, Z, r2 = fit(viscosity, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict['viscosity'] = c


VISCOSITY
Parameters: [ 6.97318785e-04 -1.47128076e-08  1.98579652e-13 -9.38312502e-19]
R2: 0.993386025713022

In [14]:
print("HEAT_CAPACITY")
c, Z_fit, Z, r2 = fit(heat_capacity, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict['heat_capacity'] = c


HEAT_CAPACITY
Parameters: [ 4.17429701e+03  4.31727962e-04  7.44863129e-10 -8.92950974e-15]
R2: 0.9998276406224097

In [15]:
print("CONDUCTIVITY")
c, Z_fit, Z, r2 = fit(conductivity, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict['conductivity'] = c


CONDUCTIVITY
Parameters: [ 6.20940155e-01  1.90022194e-06 -2.45476728e-11  1.12611502e-16]
R2: 0.9960286787923832

In [16]:
print("VAPORIZATION ENTHALPY")
c, Z_fit, Z, r2 = fit(vaporization_enthalpy, X)
plt.plot(X, Z, 'r') # plotting t, a separately 
plt.plot(X, Z_fit, 'b') # plotting t, b separately 
plt.show()
print("Parameters:",c)
print("R2:", r2)
c_dict['vaporization_enthalpy'] = c


VAPORIZATION ENTHALPY
Parameters: [ 2.42706159e+06 -4.20868251e+00  4.44518182e-05 -1.94852702e-10]
R2: 0.9990556364634905

In [17]:
import pprint
pprint.pprint(c_dict)


{'conductivity': array([ 6.20940155e-01,  1.90022194e-06, -2.45476728e-11,  1.12611502e-16]),
 'density': array([ 9.96906358e+02, -8.11519963e-04,  7.28615211e-09, -3.02080951e-14]),
 'enthalpy': array([ 1.31488947e+05,  7.28807358e+00, -7.84449119e-05,  3.44521411e-10]),
 'heat_capacity': array([ 4.17429701e+03,  4.31727962e-04,  7.44863129e-10, -8.92950974e-15]),
 'saturation_temperature': array([ 3.04528132e+02,  1.74387838e-03, -1.88042391e-08,  8.25850104e-14]),
 'vaporization_enthalpy': array([ 2.42706159e+06, -4.20868251e+00,  4.44518182e-05, -1.94852702e-10]),
 'vapour_density': array([ 4.62513512e-03,  6.47401547e-06, -9.47682219e-12,  3.32572527e-17]),
 'vapour_total_compressibility': array([ 6.47401547e-06, -1.89536444e-11,  9.97717582e-17]),
 'viscosity': array([ 6.97318785e-04, -1.47128076e-08,  1.98579652e-13, -9.38312502e-19])}

In [ ]:


In [ ]: