tnm function with numpy arrays

src/nmreval/dsc/TNMH_v0c-tau.py

@@ -29,11 +29,11 @@ def tau_k(tau_0, deltaE, temp_k, xx, T_f_km1):
 dataName = sys.argv[1]
 try:
     data = np.loadtxt(dataName, skiprows=0).T
-    print("Loading")
+    # print("Loading")
     temp = data[0]
     heat_capacity = data[1]
 except IndexError:
-    print("File not found")
+    # print("File not found")
     exit()
 
 
@@ -89,27 +89,25 @@ T_f_ns = []
 R = 8.314
 
 
-def tnmfunc(temperature, tau_g, xx, beta, deltaE):
+def tnmfunc(temperature, tau_g, x, beta, deltaE):
-    steps = len(temperature)
+    step = len(temperature)
-    tau = np.zeros(steps)
-    temp_fictive = np.zeros(steps)
 
-    # Initialize first values
+    Tf = np.empty(step)
-    temp_fictive[0] = temperature[0]
+    Tf[0] = temperature[0]
-    tau[0] = relax(temp_fictive[0], temp_fictive[0], tau_g, T_g, deltaE, xx)
 
-    deltaT = np.zeros(steps)
+    delT = np.diff(temperature)
-    deltaT[1:] = np.diff(temperature)
+    delt = np.abs(delT) * 60 / rate
-    delta_t = deltaT * 60 / rate
-    dttau = np.zeros(steps)
 
-    for i in range(1, steps):
+    tau = np.empty(step)
-        tau[i] = relax(temperature[i], temp_fictive[i-1], tau_g, T_g, deltaE, xx)
+    dttau = np.zeros(step)
-        dttau[:i] += delta_t[i-1] / tau[i]
+    temp_0 = temperature[0]
-        temp_fictive[i] = np.sum(deltaT[:i] * (1-np.exp(-(dttau[:i])**beta))) + temp_fictive[0]
-        print(i, temperature[i], temp_fictive[i-i], dttau, tau[i], delta_t[i-1], deltaT[i-1])
 
-    return temp_fictive
+    for i in range(0, step-1):
+        tau[i] = relax(temperature[i+1], Tf[i], tau_g, T_g, deltaE, x)
+        dttau[:i] += delt[i] / tau[i]
+        Tf[i+1] = np.sum(delT[:i] * (1-np.exp(-dttau[:i]**beta))) + temp_0
+
+    return Tf
 
 
 def relax(t, tf, tau_g, t_glass, ea, x):
@@ -153,7 +151,101 @@ def tnmfunc2(xdata, tau_0, xx, beta, deltaE):
         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])
+        # 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))))
@@ -188,10 +280,31 @@ def fitTNMH(x, tau_0, xx, beta, deltaE):
 
 # Use only every 20th data point to reduce fitting time
 temperaturesInterpol = np.linspace(temperatures[0], temperatures[-1], 20)
-temperaturesInterpol = np.linspace(150, 200, 5)
+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))
+plt.plot(tnmfunc2(temperaturesInterpol, *p0), label='Flo')
-plt.plot(tnmfunc(temperaturesInterpol, *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(np.diff(tnmfunc(temperaturesInterpol, *p0)), label='array')
+plt.legend()
 plt.show()
 # print(fitTNMH(temperaturesInterpol, *p0))