diff --git a/src/nmreval/dsc/TNMH_v0c-tau.py b/src/nmreval/dsc/TNMH_v0c-tau.py index 28e505f..be951fe 100644 --- a/src/nmreval/dsc/TNMH_v0c-tau.py +++ b/src/nmreval/dsc/TNMH_v0c-tau.py @@ -1,41 +1,26 @@ -# Tool-Naranayaswamy-Moynihan-Hodge model -# by Florian Pabst 2020 - - -import re -import os import sys -import time import numpy as np import matplotlib.pyplot as plt -from scipy.integrate import quad, cumulative_trapezoid +from scipy.integrate import cumulative_trapezoid from scipy.optimize import curve_fit, fsolve from scipy.stats import linregress from nmreval.data import Points -def slope(x, t, m): - return t+x*m - - -def tau_k(tau_0, deltaE, temp_k, xx, T_f_km1): - return tau_0 * np.exp(xx*deltaE / (R*temp_k) + (1-xx) * deltaE / (R*T_f_km1)) - - # Read data dataName = sys.argv[1] try: data = np.loadtxt(dataName, skiprows=0).T - # print("Loading") temp = data[0] heat_capacity = data[1] except IndexError: - # print("File not found") exit() +print(np.round(temp, decimals=3)[:30]) + def get_fictive_temperature(pts: Points, glass: tuple[float, float], liquid: tuple[float, float]): min_glass, max_glass = min(glass), max(glass) @@ -51,8 +36,8 @@ def get_fictive_temperature(pts: Points, glass: tuple[float, float], liquid: tup region.y -= glass_extrapolation liquid_regime = (min_liquid < region.x) & (region.x < max_liquid) - regress = linregress(region.x[liquid_regime], region.y[liquid_regime]) - liquid_extrapolation = regress.slope * region.x + regress.intercept + regress2 = linregress(region.x[liquid_regime], region.y[liquid_regime]) + liquid_extrapolation = regress2.slope * region.x + regress2.intercept real_area = -cumulative_trapezoid(region.y, region.x, initial=0) real_area -= real_area[-1] @@ -60,10 +45,12 @@ def get_fictive_temperature(pts: Points, glass: tuple[float, float], liquid: tup equivalent_area -= equivalent_area[-1] equivalent_area *= -1 - return region.x, region.x[np.argmin(np.abs(real_area[:, None] - equivalent_area[None, :]), axis=1)] + return region.x, np.round(region.x[np.argmin(np.abs(real_area[:, None] - equivalent_area[None, :]), axis=1)], decimals=4) curve = Points(x=temp, y=heat_capacity, name=dataName) +# plt.plot(curve.x, curve.y) +# plt.show() # First step: Calculate fictive Temperature T_f(T) from Cp data, then dT_f/dT @@ -75,21 +62,19 @@ initial_Tf = np.mean(fictiveTemp[:20]) rate = curve.value # Calculate dTf/dT -with np.errstate(all='ignore'): - dTfdT = np.gradient(fictiveTemp, temperatures) -nan_filter = ~np.isnan(dTfdT) -temperatures = np.array(temperatures)[nan_filter] -dTfdT = dTfdT[nan_filter] +# with np.errstate(all='ignore'): +# dTfdT = np.diff(fictiveTemp)/np.diff(temperatures) # np.gradient(fictiveTemp, temperatures) +# nan_filter = ~np.isnan(dTfdT) +# temperatures = 0.5 * (temperatures[:-1] + temperatures[1:])[nan_filter] +# # fictiveTemp = fictiveTemp[nan_filter] +# dTfdT = dTfdT[nan_filter] +print(fictiveTemp[:30]) -exponents = [] -taus = [] -temps = [0] -T_f_ns = [] R = 8.314 -def tnmfunc(temperature, tau_g, x, beta, deltaE): +def tnm_function(temperature, tau_g, x, beta, energy, tg): step = len(temperature) Tf = np.empty(step) @@ -103,7 +88,7 @@ def tnmfunc(temperature, tau_g, x, beta, deltaE): temp_0 = temperature[0] for i in range(0, step-1): - tau[i] = relax(temperature[i+1], Tf[i], tau_g, T_g, deltaE, x) + tau[i] = relax(temperature[i+1], Tf[i], tau_g, tg, energy, x) dttau[:i] += delt[i] / tau[i] Tf[i+1] = np.sum(delT[:i] * (1-np.exp(-dttau[:i]**beta))) + temp_0 @@ -115,147 +100,6 @@ def relax(t, tf, tau_g, t_glass, ea, x): return tau_g * np.exp((x*h / t) + ((1 - x) * h / tf) - h / t_glass) -# TNMH Code (adapted from Badrinarayanan's PhD Thesis): -def tnmfunc2(xdata, tau_0, xx, beta, deltaE): - delhR = deltaE/8.314 - T = [] - Tf = [] - t = [] - delt = [] - tau = [] - dttau = [] - delT = [] - step = len(xdata) - for k in range(0, step): - dttau.append(0) - T.append(240) - Tf.append(240) - t.append(0) - delt.append(0) - tau.append(0) - dttau.append(0) - delT.append(0) - T[1] = xdata[1] - Tf[1] = T[1] - t[1] = 0 - delt[1] = 0 - tau[1] = tau_0 * np.exp((xx*delhR / T[1]) + ((1 - xx)*delhR / Tf[1])-delhR/T_g) - dttau[1] = delt[1] / tau[1] - delT[1] = 0 - - Tfinit = T[1] - T_f_ns = [xdata[0]] - for i in range(1, step): - T[i] = xdata[i] - delT[i] = xdata[i] - xdata[i-1] - delt[i] = abs(delT[i]) / (rate/60) - t[i] = t[i-1] + delt[i] - tau[i] = tau_0 * np.exp((delhR*xx / T[i]) + ((1 - xx)*delhR / Tf[i-1])-delhR/T_g) - # print(tau[i], T[i], Tf[i-1], dttau) - for j in np.arange(2, i).reshape(-1): - dttau[j] = dttau[j] + (delt[i] / tau[i]) - Tfinit = Tfinit + (delT[j] * (1 - np.exp(- (dttau[j] ** beta)))) - Tf[i] = Tfinit - T_f_ns.append(Tfinit) - Tfinit = Tf[1] - - return T_f_ns - - -def tnmfunc3(xdata, tau_0, xx, beta, deltaE): - delhR = deltaE/8.314 - - step = len(xdata) - - dttau = [] - T = [] - Tf = [] - t = [] - delT = [] - delt = [] - tau = [] - for _ in range(step): - T.append(np.nan) - Tf.append(np.nan) - t.append(np.nan) - delT.append(np.nan) - delt.append(np.nan) - tau.append(np.nan) - dttau.append(0) - - T[0] = xdata[0] - Tf[0] = T[0] - t[0] = 0 - delT[0] = 0 - delt[0] = 0 - tau[0] = tau_0 * np.exp(xx * delhR/T[0] + (1-xx)*delhR/Tf[0] - delhR/T_g) - dttau[0] = delt[0]/tau[0] - - Tfinit = T[0] - - for i in range(1, step): - T[i] = xdata[i] - delT[i] = xdata[i] - xdata[i-1] - delt[i] = abs(delT[i]) * 60 / rate - t[i] = t[i-1] + delt[i] - tau[i] = tau_0 * np.exp(xx * delhR/T[i] + (1-xx)*delhR/Tf[i-1] - delhR/T_g) - - for j in range(1, i): - dttau[j] += delt[i]/tau[i] - Tfinit += delT[j] * (1-np.exp(-dttau[j]**beta)) - - Tf[i] = Tfinit - Tfinit = T[0] - - return Tf - - -def tnmfunc2(xdata, tau_0, xx, beta, deltaE): - delhR = deltaE/8.314 - T = [] - Tf = [] - t = [] - delt = [] - tau = [] - dttau = [] - delT = [] - step = len(xdata) - for k in range(0, step): - dttau.append(0) - T.append(240) - Tf.append(240) - t.append(0) - delt.append(0) - tau.append(0) - dttau.append(0) - delT.append(0) - T[1] = xdata[1] - Tf[1] = T[1] - t[1] = 0 - delt[1] = 0 - tau[1] = tau_0 * np.exp((xx*delhR / T[1]) + ((1 - xx)*delhR / Tf[1])-delhR/T_g) - dttau[1] = delt[1] / tau[1] - delT[1] = 0 - - Tfinit = T[1] - T_f_ns = [xdata[0]] - for i in range(1, step): - T[i] = xdata[i] - delT[i] = xdata[i] - xdata[i-1] - delt[i] = abs(delT[i]) / (rate/60) - t[i] = t[i-1] + delt[i] - tau[i] = tau_0 * np.exp((delhR*xx / T[i]) + ((1 - xx)*delhR / Tf[i-1])-delhR/T_g) - # print(tau[i], T[i], Tf[i-1], dttau) - for j in np.arange(2, i).reshape(-1): - dttau[j] = dttau[j] + (delt[i] / tau[i]) - Tfinit = Tfinit + (delT[j] * (1 - np.exp(- (dttau[j] ** beta)))) - Tf[i] = Tfinit - T_f_ns.append(Tfinit) - Tfinit = Tf[1] - - return T_f_ns - - # Use T_g = T_fictive for tau determination T_g = initial_Tf @@ -263,48 +107,26 @@ T_g = initial_Tf p0 = [10, 0.5, 0.4, 225085] -def fitTNMH(x, tau_0, xx, beta, deltaE): - dTNMHdT = [] - modelTemp = np.append(x[::-1], x[1:]) +def tnmh_fit(tg): + def wrap(x, tau_0, xx, beta, energy): + modelTemp = np.r_[x[::-1], x[1:]] - TNMH = tnmfunc2(modelTemp, tau_0, xx, beta, deltaE) - for i in range(len(modelTemp)-1): - dTNMHdT.append((TNMH[i+1]-TNMH[i]) / (modelTemp[i+1]-modelTemp[i])) - res = dTNMHdT[len(modelTemp)//2-1:] + TNMH = tnm_function(modelTemp, tau_0, xx, beta, energy, tg) + res = np.gradient(TNMH, modelTemp) - plt.plot(x, res) - plt.show() + return res[len(x)-1:] - return np.array(res) + return wrap # Use only every 20th data point to reduce fitting time -temperaturesInterpol = np.linspace(temperatures[0], temperatures[-1], 20) -temperaturesInterpol = np.linspace(150, 200, 51) - -temperaturesInterpol = np.append(temperaturesInterpol[::-1], temperaturesInterpol[1:]) - -start = time.time() -tnmfunc2(temperaturesInterpol, *p0) -print('Flo', time.time()-start) -start = time.time() -tnmfunc3(temperaturesInterpol, *p0) -print('list', time.time()-start) -start = time.time() -tnmfunc(temperaturesInterpol, *p0) -print('Array', time.time()-start) - +temperaturesInterpol = np.linspace(temperatures[0], temperatures[-1], 200) dTfdTInterpol = np.interp(temperaturesInterpol, temperatures, dTfdT) -plt.plot(tnmfunc2(temperaturesInterpol, *p0), label='Flo') -plt.plot(tnmfunc3(temperaturesInterpol, *p0), label='List') -plt.plot(tnmfunc(temperaturesInterpol, *p0), label='array') -plt.legend() -plt.show() +res = curve_fit(tnmh_fit(T_g), temperaturesInterpol, dTfdTInterpol, p0) -plt.plot(np.diff(tnmfunc2(temperaturesInterpol, *p0)), label='Flo') -plt.plot(np.diff(tnmfunc3(temperaturesInterpol, *p0)), label='List') -plt.plot(np.diff(tnmfunc(temperaturesInterpol, *p0)), label='array') + +plt.plot(temperatures, dTfdT, label='dTf/dT') +plt.plot(temperaturesInterpol, tnmh_fit(T_g)(temperaturesInterpol, *res[0]), label='fit') plt.legend() plt.show() -# print(fitTNMH(temperaturesInterpol, *p0)) diff --git a/src/nmreval/dsc/TNMH_v0c-tau2.py b/src/nmreval/dsc/TNMH_v0c-tau2.py index 3ac2e0d..6ca2216 100644 --- a/src/nmreval/dsc/TNMH_v0c-tau2.py +++ b/src/nmreval/dsc/TNMH_v0c-tau2.py @@ -26,7 +26,6 @@ except IndexError: dataOutName = os.path.splitext(str(dataName))[0] + "_Tfict+TNMH.dat" - # First step: Calculate fictive Temperature T_f(T) from Cp data, then dT_f/dT # Find start and end point for glass Cp linear fit @@ -111,7 +110,7 @@ res = curve_fit(slope, dataTemp[fitLimitLeft_Gl:fitLimitRight_Gl], dataCp[fitLim glassCpParam = res[0] # Determine liquid Cp slope -res = curve_fit(slope, dataTemp[fitLimitLeft_Lq:fitLimitRight_Lq], dataCp[fitLimitLeft_Lq:fitLimitRight_Lq]) +res = curve_fit(slope, dataTemp[fitLimitLeft_Lq:fitLimitRight_Lq], dataCp[fitLimitLeft_Lq:fitLimitRight_Lq]) liquidCpParam = res[0] xLin = np.linspace(coordsGl[0][0], coordsLq[1][0], len(dataTemp)) @@ -138,7 +137,7 @@ intStop = fitLimitRight_Lq fictiveTemp = [] temperatures = [] -for i,tempStart in enumerate(dataTemp): +for i, tempStart in enumerate(dataTemp): if fitLimitLeft_Gl < i < fitLimitRight_Lq: intStart = i @@ -172,20 +171,15 @@ for i in range(len(temperatures)-1): dTfdT.append((fictiveTemp[i+1]-fictiveTemp[i]) / (temperatures[i+1]-temperatures[i])) temperaturesCut.append(temperatures[i+1]) -""" -fig, (ax1) = P.subplots(1) +print(fictiveTemp[:30]) + +fig, ax1 = plt.subplots(1) ax1.set_xlabel('Temperature / K') ax1.set_ylabel('Fictive Temperature / K') -#ax1.set_xlim(left=coordsGl[0][0]) -#ax1.set_xlim(right=coordsLq[1][0]) -#ax1.set_ylim(bottom=coordsGl[0][1]-0.1) -#ax1.set_ylim(top=coordsLq[1][1]+0.1) -ax1.plot(temperatures,fictiveTemp, 'ro') -#ax1.plot(temperatures,(T_f_0+temperatures*0), 'b-') -ax1.axhline(y=T_f_0[0], color='b', linestyle='-') -#ax1.plot(temperatures[:-1],dTfdT, 'ro') -P.show() -""" +# ax1.plot(temperatures, fictiveTemp, 'ro') +ax1.plot(temperaturesCut, dTfdT) +# ax1.axhline(y=T_f_0[0], color='b', linestyle='-') +plt.show() # Second step: Calcualte and Fit TNMH-model try: @@ -207,6 +201,7 @@ T_f_ns = [] markovs = [coordsLq[1][0]] R = 8.31 + # TNMH Code (adapted from Badrinarayanan's PhD Thesis): def tnmfunc2(xdata, tau_0, xx, beta, deltaE): delhR = deltaE/8.314 @@ -250,7 +245,6 @@ def tnmfunc2(xdata, tau_0, xx, beta, deltaE): T_f_ns.append(Tfinit) Tfinit = Tf[1] - print(T_f_ns[:10]) return T_f_ns @@ -297,25 +291,15 @@ tau_0, xx, beta, deltaE = p0 fitResult = fitTNMH(p0, temperaturesInterpol) -""" -fig, (ax1) = P.subplots(1) + +fig, ax1 = plt.subplots() ax1.set_xlabel('Temperature / K') ax1.set_ylabel('Fictive Temperature / K') -#ax1.set_xlim(left=coordsGl[0][0]) -#ax1.set_xlim(right=coordsLq[1][0]) -#ax1.set_ylim(bottom=coordsGl[0][1]-0.1) -#ax1.set_ylim(top=coordsLq[1][1]+0.1) -#ax1.plot(temperatures,fictiveTemp, 'bo') -#ax1.plot(temperatures,T_f_ns, 'r-') -ax1.plot(temperaturesCut,dTfdT, 'bo') -#ax1.plot(temperatures[0:5454],dTNMHdT[0:5454], 'bo') -#ax1.plot(temperatures[5456:],dTNMHdT[5455:], 'ro') -#ax1.plot(temperaturesModel[0:99],dTNMHdT[0:99], 'b-') -#ax1.plot(temperaturesModel[100:199],dTNMHdT[100:199], 'r-') -ax1.plot(temperaturesInterpol,dTfdTInterpol, 'yo') -ax1.plot(temperaturesInterpol,fitResult, 'r-') -P.show() -""" +ax1.plot(temperaturesCut, dTfdT, 'bo') +ax1.plot(temperaturesInterpol, dTfdTInterpol, 'yo') +ax1.plot(temperaturesInterpol, fitResult, 'r-') +plt.show() + dataOutName = os.path.splitext(str(dataName))[0] + "_dTfdT_tau.fit" dataOut = np.column_stack((temperaturesInterpol, dTfdTInterpol, fitResult))