In [2]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import os,sys,inspect
currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
parentdir = os.path.dirname(currentdir)
sys.path.insert(0,parentdir)
from water_properties import density, enthalpy, heat_capacity, conductivity, viscosity
In [3]:
def simplified_model_for_fit(x, a00, a01, a02, a10, a11, a12, a20, a21, a22):
c = [[a00, a01, a02],[a10, a11, a12],[a20, a21, a22]]
T, P = x
return simplified_model(T, P, c)
In [4]:
def simplified_model(T, P, c):
return np.polynomial.polynomial.polyval2d(T, P, c)
In [5]:
def calculate_output(func,X1,X2,*args):
Z = np.zeros_like(X1)
for i in range(X1.shape[0]):
for j in range(X2.shape[1]):
Z[i,j] = func(X1[i,j],X2[i,j],*args)[0]
return Z
In [15]:
limits = [278., 353., 10000., 1000000.]
npoints = 100
T = np.linspace(limits[0], limits[1], npoints)
P = np.linspace(limits[2], limits[3], npoints)
X1, X2 = np.meshgrid(T, P)
In [16]:
def fit(model, X1, X2):
Z = calculate_output(model,X1,X2)
size = X1.shape
x1_1d = X1.reshape((1, np.prod(size)))
x2_1d = X2.reshape((1, np.prod(size)))
z_1d = Z.reshape((1, np.prod(size)))
xdata = np.vstack((x1_1d, x2_1d))
ydata = z_1d[0]
popt, pcov = curve_fit(simplified_model_for_fit, xdata, ydata)
c = popt.reshape((3,3))
Z_fit = simplified_model(X1,X2,c)
z_fit_1d = Z_fit.reshape((1, np.prod(size)))
# residual sum of squares
ss_res = np.sum((z_1d[0] - z_fit_1d[0]) ** 2)
# total sum of squares
ss_tot = np.sum((z_1d[0] - np.mean(z_1d[0])) ** 2)
# r-squared
r2 = 1 - (ss_res / ss_tot)
return (c, Z_fit, Z, r2)
In [17]:
def plot_result(name, X1, X2, Z, Z_fit):
plt.subplot(1, 2, 1)
plt.title("Real {0}".format(name))
plt.pcolormesh(X1, X2, Z)
plt.axis(limits)
plt.colorbar()
plt.subplot(1, 2, 2)
plt.title("Fitted {0}".format(name))
plt.pcolormesh(X1, X2, Z_fit)
plt.axis(limits)
plt.colorbar()
plt.show()
In [18]:
print("DENSITY")
c, Z_fit, Z, r2 = fit(density, X1, X2)
plot_result("", X1, X2, Z, Z_fit)
print("Parameters:",c)
print("R2:", r2)
c_dict = dict()
c_dict['density'] = c
In [19]:
print("ENTHALPY")
c, Z_fit, Z, r2 = fit(enthalpy, X1, X2)
plot_result("", X1, X2, Z, Z_fit)
print("Parameters:",c)
print("R2:", r2)
c_dict['enthalpy'] = c
In [20]:
print("VISCOSITY")
c, Z_fit, Z, r2 = fit(viscosity, X1, X2)
plot_result("", X1, X2, Z, Z_fit)
print("Parameters:",c)
print("R2:", r2)
c_dict['viscosity'] = c
In [21]:
print("HEAT CAPACITY")
c, Z_fit, Z, r2 = fit(heat_capacity, X1, X2)
plot_result("", X1, X2, Z, Z_fit)
print("Parameters:",c)
print("R2:", r2)
c_dict['heat_capacity'] = c
In [22]:
print("CONDUCTIVITY")
c, Z_fit, Z, r2 = fit(conductivity, X1, X2)
plot_result("", X1, X2, Z, Z_fit)
print("Parameters:",c)
print("R2:", r2)
c_dict['conductivity'] = c
In [23]:
import pprint
pprint.pprint(c_dict)
In [ ]:
In [ ]: