some stuff done

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])
+    regress2 = linregress(region.x[liquid_regime], region.y[liquid_regime])
-    liquid_extrapolation = regress.slope * region.x + regress.intercept
+    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'):
+# with np.errstate(all='ignore'):
-    dTfdT = np.gradient(fictiveTemp, temperatures)
+#     dTfdT = np.diff(fictiveTemp)/np.diff(temperatures)  #  np.gradient(fictiveTemp, temperatures)
-nan_filter = ~np.isnan(dTfdT)
+# nan_filter = ~np.isnan(dTfdT)
-temperatures = np.array(temperatures)[nan_filter]
+# temperatures = 0.5 * (temperatures[:-1] + temperatures[1:])[nan_filter]
-dTfdT = dTfdT[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):
+def tnmh_fit(tg):
-    dTNMHdT = []
+    def wrap(x, tau_0, xx, beta, energy):
-    modelTemp = np.append(x[::-1], x[1:])
+        modelTemp = np.r_[x[::-1], x[1:]]
 
-    TNMH = tnmfunc2(modelTemp, tau_0, xx, beta, deltaE)
+        TNMH = tnm_function(modelTemp, tau_0, xx, beta, energy, tg)
-    for i in range(len(modelTemp)-1):
+        res = np.gradient(TNMH, modelTemp)
-        dTNMHdT.append((TNMH[i+1]-TNMH[i]) / (modelTemp[i+1]-modelTemp[i]))
-    res = dTNMHdT[len(modelTemp)//2-1:]
 
-    plt.plot(x, res)
+        return res[len(x)-1:]
-    plt.show()
 
-    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(temperatures[0], temperatures[-1], 200)
-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)
-
 dTfdTInterpol = np.interp(temperaturesInterpol, temperatures, dTfdT)
 
-plt.plot(tnmfunc2(temperaturesInterpol, *p0), label='Flo')
+res = curve_fit(tnmh_fit(T_g), temperaturesInterpol, dTfdTInterpol, p0)
-plt.plot(tnmfunc3(temperaturesInterpol, *p0), label='List')
-plt.plot(tnmfunc(temperaturesInterpol, *p0), label='array')
-plt.legend()
-plt.show()
 
-plt.plot(np.diff(tnmfunc2(temperaturesInterpol, *p0)), label='Flo')
+
-plt.plot(np.diff(tnmfunc3(temperaturesInterpol, *p0)), label='List')
+plt.plot(temperatures, dTfdT, label='dTf/dT')
-plt.plot(np.diff(tnmfunc(temperaturesInterpol, *p0)), label='array')
+plt.plot(temperaturesInterpol, tnmh_fit(T_g)(temperaturesInterpol, *res[0]), label='fit')
 plt.legend()
 plt.show()
-# print(fitTNMH(temperaturesInterpol, *p0))

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])
 
-"""        
+print(fictiveTemp[:30])
-fig, (ax1) = P.subplots(1)
+
+fig, ax1 = plt.subplots(1)
 ax1.set_xlabel('Temperature / K')
 ax1.set_ylabel('Fictive Temperature / K')
-#ax1.set_xlim(left=coordsGl[0][0])
+# ax1.plot(temperatures, fictiveTemp, 'ro')
-#ax1.set_xlim(right=coordsLq[1][0])
+ax1.plot(temperaturesCut, dTfdT)
-#ax1.set_ylim(bottom=coordsGl[0][1]-0.1)
+# ax1.axhline(y=T_f_0[0], color='b', linestyle='-')
-#ax1.set_ylim(top=coordsLq[1][1]+0.1)
+plt.show()
-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()
-"""
 
 # 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.plot(temperaturesCut, dTfdT, 'bo')
-#ax1.set_xlim(right=coordsLq[1][0])
+ax1.plot(temperaturesInterpol, dTfdTInterpol, 'yo')
-#ax1.set_ylim(bottom=coordsGl[0][1]-0.1)
+ax1.plot(temperaturesInterpol, fitResult, 'r-')
-#ax1.set_ylim(top=coordsLq[1][1]+0.1)
+plt.show()
-#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()
-"""
 dataOutName = os.path.splitext(str(dataName))[0] + "_dTfdT_tau.fit"
 dataOut = np.column_stack((temperaturesInterpol, dTfdTInterpol, fitResult))