add tnmh stuff
This commit is contained in:
parent
15c959fd71
commit
11ed3c16eb
40067
src/nmreval/dsc/D84_10.0K-min_h.dat
Normal file
40067
src/nmreval/dsc/D84_10.0K-min_h.dat
Normal file
File diff suppressed because it is too large
Load Diff
202
src/nmreval/dsc/D84_10.0K-min_h_dTfdT_tau.fit
Normal file
202
src/nmreval/dsc/D84_10.0K-min_h_dTfdT_tau.fit
Normal file
@ -0,0 +1,202 @@
|
||||
# tau_TNMH=77.93669386254352 beta=0.6334823578561849 x=0.4107577725462226 deltaE=205554.80338165778
|
||||
# Temp dTfdT_data dTfdT_fit
|
||||
1.50007e+02 -1.54659e-03 5.04753e-03
|
||||
1.50258e+02 -1.35164e-03 -5.81695e-03
|
||||
1.50509e+02 -1.15321e-03 -5.24594e-03
|
||||
1.50760e+02 1.33378e-03 -5.19023e-03
|
||||
1.51011e+02 -4.19078e-03 -5.23034e-03
|
||||
1.51262e+02 -3.69741e-03 -5.32106e-03
|
||||
1.51514e+02 -1.17892e-03 -5.44569e-03
|
||||
1.51765e+02 -4.26847e-04 -5.59617e-03
|
||||
1.52016e+02 1.34823e-03 -5.76806e-03
|
||||
1.52267e+02 -8.53837e-04 -5.95864e-03
|
||||
1.52518e+02 -3.32957e-03 -6.16618e-03
|
||||
1.52769e+02 -3.18265e-03 -6.38951e-03
|
||||
1.53021e+02 -9.27941e-04 -6.62778e-03
|
||||
1.53272e+02 2.81209e-03 -6.88034e-03
|
||||
1.53523e+02 2.53703e-03 -7.14668e-03
|
||||
1.53774e+02 2.79565e-03 -7.42632e-03
|
||||
1.54025e+02 3.73709e-03 -7.71879e-03
|
||||
1.54277e+02 -4.74193e-04 -8.02362e-03
|
||||
1.54528e+02 5.53513e-03 -8.34027e-03
|
||||
1.54779e+02 5.40911e-03 -8.66813e-03
|
||||
1.55030e+02 4.61145e-03 -9.00648e-03
|
||||
1.55281e+02 3.12182e-03 -9.35449e-03
|
||||
1.55532e+02 6.75992e-03 -9.71118e-03
|
||||
1.55784e+02 3.14577e-03 -1.00754e-02
|
||||
1.56035e+02 2.59412e-03 -1.04458e-02
|
||||
1.56286e+02 -1.54718e-03 -1.08208e-02
|
||||
1.56537e+02 1.17037e-04 -1.11987e-02
|
||||
1.56788e+02 -2.03614e-03 -1.15774e-02
|
||||
1.57039e+02 1.33113e-03 -1.19545e-02
|
||||
1.57291e+02 -5.96949e-03 -1.23274e-02
|
||||
1.57542e+02 -6.48426e-03 -1.26931e-02
|
||||
1.57793e+02 -1.59997e-03 -1.30481e-02
|
||||
1.58044e+02 -3.72898e-04 -1.33886e-02
|
||||
1.58295e+02 1.30990e-03 -1.37104e-02
|
||||
1.58547e+02 4.08780e-03 -1.40088e-02
|
||||
1.58798e+02 -7.05642e-03 -1.42784e-02
|
||||
1.59049e+02 -2.00584e-03 -1.45134e-02
|
||||
1.59300e+02 -1.92937e-03 -1.47073e-02
|
||||
1.59551e+02 -6.58517e-04 -1.48531e-02
|
||||
1.59802e+02 7.32033e-04 -1.49429e-02
|
||||
1.60054e+02 2.36722e-03 -1.49682e-02
|
||||
1.60305e+02 -1.17443e-03 -1.49196e-02
|
||||
1.60556e+02 4.59116e-04 -1.47868e-02
|
||||
1.60807e+02 8.17255e-03 -1.45588e-02
|
||||
1.61058e+02 8.48647e-03 -1.42234e-02
|
||||
1.61310e+02 7.96560e-03 -1.37674e-02
|
||||
1.61561e+02 1.29202e-02 -1.31766e-02
|
||||
1.61812e+02 4.84645e-03 -1.24355e-02
|
||||
1.62063e+02 1.24679e-02 -1.15274e-02
|
||||
1.62314e+02 1.08789e-02 -1.04342e-02
|
||||
1.62565e+02 9.31453e-03 -9.13630e-03
|
||||
1.62817e+02 1.10428e-02 -7.61265e-03
|
||||
1.63068e+02 1.81300e-02 -5.84047e-03
|
||||
1.63319e+02 9.70774e-03 -3.79521e-03
|
||||
1.63570e+02 1.20566e-02 -1.45035e-03
|
||||
1.63821e+02 1.49584e-02 1.22268e-03
|
||||
1.64072e+02 1.72350e-02 4.25478e-03
|
||||
1.64324e+02 2.03471e-02 7.67928e-03
|
||||
1.64575e+02 3.74633e-02 1.15322e-02
|
||||
1.64826e+02 2.20722e-02 1.58527e-02
|
||||
1.65077e+02 2.90447e-02 2.06828e-02
|
||||
1.65328e+02 3.56891e-02 2.60686e-02
|
||||
1.65580e+02 4.02175e-02 3.20598e-02
|
||||
1.65831e+02 4.99873e-02 3.87108e-02
|
||||
1.66082e+02 6.52373e-02 4.60809e-02
|
||||
1.66333e+02 5.16469e-02 5.42346e-02
|
||||
1.66584e+02 6.50874e-02 6.32432e-02
|
||||
1.66835e+02 7.59582e-02 7.31846e-02
|
||||
1.67087e+02 8.82617e-02 8.41446e-02
|
||||
1.67338e+02 9.75569e-02 9.62183e-02
|
||||
1.67589e+02 1.08193e-01 1.09511e-01
|
||||
1.67840e+02 1.09727e-01 1.24138e-01
|
||||
1.68091e+02 1.30829e-01 1.40230e-01
|
||||
1.68342e+02 1.55907e-01 1.57932e-01
|
||||
1.68594e+02 1.71526e-01 1.77404e-01
|
||||
1.68845e+02 1.87568e-01 1.98829e-01
|
||||
1.69096e+02 2.05094e-01 2.22407e-01
|
||||
1.69347e+02 2.27806e-01 2.48367e-01
|
||||
1.69598e+02 2.58337e-01 2.76962e-01
|
||||
1.69850e+02 2.88510e-01 3.08478e-01
|
||||
1.70101e+02 3.20528e-01 3.43232e-01
|
||||
1.70352e+02 3.57173e-01 3.81577e-01
|
||||
1.70603e+02 4.02201e-01 4.23903e-01
|
||||
1.70854e+02 4.47497e-01 4.70630e-01
|
||||
1.71105e+02 5.06370e-01 5.22210e-01
|
||||
1.71357e+02 5.66666e-01 5.79104e-01
|
||||
1.71608e+02 6.33707e-01 6.41763e-01
|
||||
1.71859e+02 7.08395e-01 7.10580e-01
|
||||
1.72110e+02 7.86941e-01 7.85817e-01
|
||||
1.72361e+02 8.80202e-01 8.67486e-01
|
||||
1.72612e+02 9.72949e-01 9.55178e-01
|
||||
1.72864e+02 1.07283e+00 1.04780e+00
|
||||
1.73115e+02 1.17098e+00 1.14323e+00
|
||||
1.73366e+02 1.26581e+00 1.23794e+00
|
||||
1.73617e+02 1.36284e+00 1.32657e+00
|
||||
1.73868e+02 1.40164e+00 1.40186e+00
|
||||
1.74120e+02 1.45067e+00 1.45520e+00
|
||||
1.74371e+02 1.46476e+00 1.47830e+00
|
||||
1.74622e+02 1.44706e+00 1.46612e+00
|
||||
1.74873e+02 1.40520e+00 1.42006e+00
|
||||
1.75124e+02 1.35811e+00 1.34934e+00
|
||||
1.75375e+02 1.26322e+00 1.26857e+00
|
||||
1.75627e+02 1.20703e+00 1.19227e+00
|
||||
1.75878e+02 1.15871e+00 1.12985e+00
|
||||
1.76129e+02 1.12062e+00 1.08413e+00
|
||||
1.76380e+02 1.09118e+00 1.05315e+00
|
||||
1.76631e+02 1.07111e+00 1.03314e+00
|
||||
1.76882e+02 1.05191e+00 1.02044e+00
|
||||
1.77134e+02 1.04878e+00 1.01241e+00
|
||||
1.77385e+02 1.04027e+00 1.00731e+00
|
||||
1.77636e+02 1.03275e+00 1.00413e+00
|
||||
1.77887e+02 1.02782e+00 1.00219e+00
|
||||
1.78138e+02 1.02566e+00 1.00108e+00
|
||||
1.78390e+02 1.02147e+00 1.00049e+00
|
||||
1.78641e+02 1.01938e+00 1.00020e+00
|
||||
1.78892e+02 1.02000e+00 1.00007e+00
|
||||
1.79143e+02 1.01835e+00 1.00002e+00
|
||||
1.79394e+02 1.01657e+00 1.00001e+00
|
||||
1.79645e+02 1.01941e+00 1.00000e+00
|
||||
1.79897e+02 1.00523e+00 1.00000e+00
|
||||
1.80148e+02 1.00689e+00 1.00000e+00
|
||||
1.80399e+02 1.00730e+00 1.00000e+00
|
||||
1.80650e+02 1.00768e+00 1.00000e+00
|
||||
1.80901e+02 1.00746e+00 1.00000e+00
|
||||
1.81152e+02 1.00606e+00 1.00000e+00
|
||||
1.81404e+02 1.00202e+00 1.00000e+00
|
||||
1.81655e+02 9.93205e-01 1.00000e+00
|
||||
1.81906e+02 9.98851e-01 1.00000e+00
|
||||
1.82157e+02 1.00478e+00 1.00000e+00
|
||||
1.82408e+02 1.00147e+00 1.00000e+00
|
||||
1.82660e+02 1.00747e+00 1.00000e+00
|
||||
1.82911e+02 9.83522e-01 1.00000e+00
|
||||
1.83162e+02 9.90690e-01 1.00000e+00
|
||||
1.83413e+02 9.94175e-01 1.00000e+00
|
||||
1.83664e+02 9.95883e-01 1.00000e+00
|
||||
1.83915e+02 9.92199e-01 1.00000e+00
|
||||
1.84167e+02 9.90755e-01 1.00000e+00
|
||||
1.84418e+02 9.86480e-01 1.00000e+00
|
||||
1.84669e+02 9.93964e-01 1.00000e+00
|
||||
1.84920e+02 9.94745e-01 1.00000e+00
|
||||
1.85171e+02 9.95895e-01 1.00000e+00
|
||||
1.85422e+02 9.97575e-01 1.00000e+00
|
||||
1.85674e+02 9.99326e-01 1.00000e+00
|
||||
1.85925e+02 9.97579e-01 1.00000e+00
|
||||
1.86176e+02 1.00187e+00 1.00000e+00
|
||||
1.86427e+02 9.99647e-01 1.00000e+00
|
||||
1.86678e+02 1.00026e+00 1.00000e+00
|
||||
1.86930e+02 1.00307e+00 1.00000e+00
|
||||
1.87181e+02 1.01192e+00 1.00000e+00
|
||||
1.87432e+02 9.96247e-01 1.00000e+00
|
||||
1.87683e+02 9.98530e-01 1.00000e+00
|
||||
1.87934e+02 1.00165e+00 1.00000e+00
|
||||
1.88185e+02 1.00562e+00 1.00000e+00
|
||||
1.88437e+02 1.01337e+00 1.00000e+00
|
||||
1.88688e+02 1.01309e+00 1.00000e+00
|
||||
1.88939e+02 9.95564e-01 1.00000e+00
|
||||
1.89190e+02 1.00160e+00 1.00000e+00
|
||||
1.89441e+02 1.00596e+00 1.00000e+00
|
||||
1.89692e+02 1.01222e+00 1.00000e+00
|
||||
1.89944e+02 1.00825e+00 1.00000e+00
|
||||
1.90195e+02 1.00603e+00 1.00000e+00
|
||||
1.90446e+02 1.00770e+00 1.00000e+00
|
||||
1.90697e+02 9.89746e-01 1.00000e+00
|
||||
1.90948e+02 1.00511e+00 1.00000e+00
|
||||
1.91200e+02 1.00621e+00 1.00000e+00
|
||||
1.91451e+02 1.00143e+00 1.00000e+00
|
||||
1.91702e+02 1.00224e+00 1.00000e+00
|
||||
1.91953e+02 1.03274e+00 1.00000e+00
|
||||
1.92204e+02 9.87474e-01 1.00000e+00
|
||||
1.92455e+02 9.93631e-01 1.00000e+00
|
||||
1.92707e+02 9.96782e-01 1.00000e+00
|
||||
1.92958e+02 9.99014e-01 1.00000e+00
|
||||
1.93209e+02 1.00251e+00 1.00000e+00
|
||||
1.93460e+02 1.01135e+00 1.00000e+00
|
||||
1.93711e+02 9.92486e-01 1.00000e+00
|
||||
1.93962e+02 9.94748e-01 1.00000e+00
|
||||
1.94214e+02 9.96729e-01 1.00000e+00
|
||||
1.94465e+02 9.99500e-01 1.00000e+00
|
||||
1.94716e+02 1.00314e+00 1.00000e+00
|
||||
1.94967e+02 1.00767e+00 1.00000e+00
|
||||
1.95218e+02 9.92280e-01 1.00000e+00
|
||||
1.95470e+02 9.93172e-01 1.00000e+00
|
||||
1.95721e+02 9.99065e-01 1.00000e+00
|
||||
1.95972e+02 1.00089e+00 1.00000e+00
|
||||
1.96223e+02 1.00147e+00 1.00000e+00
|
||||
1.96474e+02 1.00277e+00 1.00000e+00
|
||||
1.96725e+02 9.93517e-01 1.00000e+00
|
||||
1.96977e+02 9.89833e-01 1.00000e+00
|
||||
1.97228e+02 1.00075e+00 1.00000e+00
|
||||
1.97479e+02 1.00224e+00 1.00000e+00
|
||||
1.97730e+02 9.97354e-01 1.00000e+00
|
||||
1.97981e+02 9.96498e-01 1.00000e+00
|
||||
1.98232e+02 1.01632e+00 1.00000e+00
|
||||
1.98484e+02 9.86819e-01 1.00000e+00
|
||||
1.98735e+02 1.00374e+00 1.00000e+00
|
||||
1.98986e+02 1.00217e+00 1.00000e+00
|
||||
1.99237e+02 1.00089e+00 1.00000e+00
|
||||
1.99488e+02 1.00215e+00 1.00000e+00
|
||||
1.99740e+02 1.01197e+00 1.00000e+00
|
||||
1.99991e+02 9.92503e-01 1.00000e+00
|
7
src/nmreval/dsc/D84_DSC_tau.param
Normal file
7
src/nmreval/dsc/D84_DSC_tau.param
Normal file
@ -0,0 +1,7 @@
|
||||
|
||||
#Heating Rate tau_TNMH beta x deltaE Tg fictive
|
||||
1.0000e+01 9.0613e+01 6.2053e-01 3.9548e-01 2.0820e+05 1.6911e+02
|
||||
#Heating Rate tau_TNMH beta x deltaE Tg fictive
|
||||
1.0000e+01 3.5330e+01 6.7660e-01 5.2984e-01 1.7362e+05 1.6934e+02
|
||||
#Heating Rate tau_TNMH beta x deltaE Tg fictive
|
||||
1.0000e+01 7.7937e+01 6.3348e-01 4.1076e-01 2.0555e+05 1.6925e+02
|
197
src/nmreval/dsc/TNMH_v0c-tau.py
Normal file
197
src/nmreval/dsc/TNMH_v0c-tau.py
Normal file
@ -0,0 +1,197 @@
|
||||
# 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.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()
|
||||
|
||||
|
||||
def get_fictive_temperature(pts: Points, glass: tuple[float, float], liquid: tuple[float, float]):
|
||||
min_glass, max_glass = min(glass), max(glass)
|
||||
min_liquid, max_liquid = min(liquid), max(liquid)
|
||||
|
||||
region = pts.copy()
|
||||
region.cut(min_glass, max_liquid)
|
||||
|
||||
glass_regime = (min_glass < region.x) & (region.x < max_glass)
|
||||
regress = linregress(region.x[glass_regime], region.y[glass_regime])
|
||||
glass_extrapolation = regress.slope * region.x + regress.intercept
|
||||
|
||||
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
|
||||
|
||||
real_area = -cumulative_trapezoid(region.y, region.x, initial=0)
|
||||
real_area -= real_area[-1]
|
||||
equivalent_area = cumulative_trapezoid(liquid_extrapolation, region.x, initial=0)
|
||||
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)]
|
||||
|
||||
|
||||
curve = Points(x=temp, y=heat_capacity, name=dataName)
|
||||
|
||||
# First step: Calculate fictive Temperature T_f(T) from Cp data, then dT_f/dT
|
||||
|
||||
temperatures, fictiveTemp = get_fictive_temperature(curve, glass=(150, 160), liquid=(190, 200))
|
||||
|
||||
# Determine limiting T_f, i.e. T_f at lowest temperature
|
||||
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]
|
||||
|
||||
|
||||
exponents = []
|
||||
taus = []
|
||||
temps = [0]
|
||||
T_f_ns = []
|
||||
R = 8.314
|
||||
|
||||
|
||||
def tnmfunc(temperature, tau_g, xx, beta, deltaE):
|
||||
steps = len(temperature)
|
||||
tau = np.zeros(steps)
|
||||
temp_fictive = np.zeros(steps)
|
||||
|
||||
# Initialize first values
|
||||
temp_fictive[0] = temperature[0]
|
||||
tau[0] = relax(temp_fictive[0], temp_fictive[0], tau_g, T_g, deltaE, xx)
|
||||
|
||||
deltaT = np.zeros(steps)
|
||||
deltaT[1:] = np.diff(temperature)
|
||||
delta_t = deltaT * 60 / rate
|
||||
dttau = np.zeros(steps)
|
||||
|
||||
for i in range(1, steps):
|
||||
tau[i] = relax(temperature[i], temp_fictive[i-1], tau_g, T_g, deltaE, xx)
|
||||
dttau[:i] += delta_t[i-1] / tau[i]
|
||||
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
|
||||
|
||||
|
||||
def relax(t, tf, tau_g, t_glass, ea, x):
|
||||
h = ea/R
|
||||
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])
|
||||
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
|
||||
|
||||
# Start values
|
||||
p0 = [10, 0.5, 0.4, 225085]
|
||||
|
||||
|
||||
def fitTNMH(x, tau_0, xx, beta, deltaE):
|
||||
dTNMHdT = []
|
||||
modelTemp = np.append(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:]
|
||||
|
||||
plt.plot(x, res)
|
||||
plt.show()
|
||||
|
||||
return np.array(res)
|
||||
|
||||
|
||||
# Use only every 20th data point to reduce fitting time
|
||||
temperaturesInterpol = np.linspace(temperatures[0], temperatures[-1], 20)
|
||||
temperaturesInterpol = np.linspace(150, 200, 5)
|
||||
dTfdTInterpol = np.interp(temperaturesInterpol, temperatures, dTfdT)
|
||||
|
||||
plt.plot(tnmfunc2(temperaturesInterpol, *p0))
|
||||
plt.plot(tnmfunc(temperaturesInterpol, *p0))
|
||||
plt.show()
|
||||
# print(fitTNMH(temperaturesInterpol, *p0))
|
333
src/nmreval/dsc/TNMH_v0c-tau2.py
Normal file
333
src/nmreval/dsc/TNMH_v0c-tau2.py
Normal file
@ -0,0 +1,333 @@
|
||||
# Tool-Naranayaswami-Moynihan-Hodge model
|
||||
# by Florian Pabst 2020
|
||||
|
||||
|
||||
import time
|
||||
import re
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import scipy.odr as odr
|
||||
from scipy.integrate import quad
|
||||
from scipy.optimize import curve_fit, fsolve
|
||||
|
||||
# Read data
|
||||
dataName = sys.argv[1]
|
||||
try:
|
||||
data = np.loadtxt(dataName,skiprows=0)
|
||||
print("Loading")
|
||||
dataTemp = data[:, 0]
|
||||
dataCp = data[:, 1]
|
||||
except IndexError:
|
||||
print("File not found")
|
||||
exit()
|
||||
|
||||
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
|
||||
coordsGl = [(150, -0.126), (160, -0.09)]
|
||||
"""
|
||||
print("Please click on the left limit for glass Cp linear fit")
|
||||
fig, (ax1) = P.subplots(1)
|
||||
ax1.set_xlabel('Temperature / K')
|
||||
ax1.set_xlim(left=140)
|
||||
ax1.set_xlim(right=200)
|
||||
ax1.plot(dataTemp,dataCp, 'ro')
|
||||
#ax1.legend()
|
||||
def onclick(event):
|
||||
ix = event.xdata
|
||||
iy = event.ydata
|
||||
print(ix)
|
||||
|
||||
global coordsGl
|
||||
coordsGl.append((ix, iy))
|
||||
|
||||
if len(coordsGl) == 1:
|
||||
print('right limit')
|
||||
|
||||
if len(coordsGl) == 2:
|
||||
fig.canvas.mpl_disconnect(cid)
|
||||
|
||||
return coordsGl
|
||||
cid = fig.canvas.mpl_connect('button_press_event', onclick)
|
||||
P.show()
|
||||
"""
|
||||
print(coordsGl)
|
||||
|
||||
for i in range(0,len(dataTemp)):
|
||||
if dataTemp[i] < coordsGl[0][0]:
|
||||
fitLimitLeft_Gl = i
|
||||
if dataTemp[i] < coordsGl[1][0]:
|
||||
fitLimitRight_Gl = i
|
||||
|
||||
# Find start and end point for liquid Cp linear fit
|
||||
coordsLq = [(190, 2.07), (200, 1.773)]
|
||||
"""
|
||||
print("Please click on the left limit for liquid Cp linear fit")
|
||||
fig, (ax1) = P.subplots(1)
|
||||
ax1.set_xlabel('Temperature / K')
|
||||
ax1.set_xlim(left=180)
|
||||
ax1.set_xlim(right=250)
|
||||
ax1.plot(dataTemp,dataCp, 'ro')
|
||||
#ax1.legend()
|
||||
def onclick(event):
|
||||
ix = event.xdata
|
||||
iy = event.ydata
|
||||
print(ix)
|
||||
|
||||
global coordsLq
|
||||
coordsLq.append((ix, iy))
|
||||
|
||||
if len(coordsLq) == 1:
|
||||
print('right limit')
|
||||
|
||||
if len(coordsLq) == 2:
|
||||
fig.canvas.mpl_disconnect(cid)
|
||||
|
||||
return coordsLq
|
||||
cid = fig.canvas.mpl_connect('button_press_event', onclick)
|
||||
P.show()
|
||||
#print(coordsLq)
|
||||
"""
|
||||
|
||||
for i in range(0, len(dataTemp)):
|
||||
if dataTemp[i] < coordsLq[0][0]:
|
||||
fitLimitLeft_Lq = i
|
||||
if dataTemp[i] < coordsLq[1][0]:
|
||||
fitLimitRight_Lq = i
|
||||
|
||||
|
||||
def slope(x, t, m):
|
||||
return t+x*m
|
||||
|
||||
|
||||
# Determine glass Cp slope
|
||||
res = curve_fit(slope, dataTemp[fitLimitLeft_Gl:fitLimitRight_Gl], dataCp[fitLimitLeft_Gl:fitLimitRight_Gl])
|
||||
glassCpParam = res[0]
|
||||
|
||||
# Determine liquid Cp slope
|
||||
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))
|
||||
|
||||
CpInterpol = np.interp(xLin, dataTemp, dataCp)
|
||||
"""
|
||||
fig, (ax1) = P.subplots(1)
|
||||
ax1.set_xlabel('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(dataTemp,dataCp, 'ro')
|
||||
ax1.plot(dataTemp,slope(dataTemp,glassCpParam[0],glassCpParam[1]), 'b-')
|
||||
ax1.plot(dataTemp,slope(dataTemp,liquidCpParam[0],liquidCpParam[1]), 'b-')
|
||||
ax1.plot(xLin,CpInterpol, 'g-')
|
||||
P.show()
|
||||
"""
|
||||
|
||||
|
||||
# Calculate Fictive Temperature for all temperatures in range
|
||||
intStart = fitLimitLeft_Gl
|
||||
intStop = fitLimitRight_Lq
|
||||
|
||||
fictiveTemp = []
|
||||
temperatures = []
|
||||
for i,tempStart in enumerate(dataTemp):
|
||||
if fitLimitLeft_Gl < i < fitLimitRight_Lq:
|
||||
intStart = i
|
||||
|
||||
integration2 = np.trapz((dataCp[intStart:intStop]-slope(dataTemp[intStart:intStop], glassCpParam[0], glassCpParam[1])), dataTemp[intStart:intStop])
|
||||
|
||||
def integrand1(temp):
|
||||
return slope(temp, liquidCpParam[0], liquidCpParam[1])-slope(temp, glassCpParam[0], glassCpParam[1])
|
||||
|
||||
def fictiveTempEq(x):
|
||||
return quad(integrand1, x, coordsLq[1][0])[0] - integration2
|
||||
|
||||
ficTemp = fsolve(fictiveTempEq, np.array([160]))
|
||||
fictiveTemp.append(ficTemp[0])
|
||||
temperatures.append(tempStart)
|
||||
|
||||
|
||||
# Determine limiting T_f, i.e. T_f at lowest temperatures
|
||||
def constant(x, t):
|
||||
return t
|
||||
|
||||
|
||||
res = curve_fit(constant, temperatures[0:20], fictiveTemp[0:20])
|
||||
T_f_0 = res[0]
|
||||
print(T_f_0[0])
|
||||
|
||||
# Calculate dTf/dT
|
||||
dTfdT = []
|
||||
temperaturesCut = []
|
||||
for i in range(len(temperatures)-1):
|
||||
if (temperatures[i+1]-temperatures[i]) != 0:
|
||||
dTfdT.append((fictiveTemp[i+1]-fictiveTemp[i]) / (temperatures[i+1]-temperatures[i]))
|
||||
temperaturesCut.append(temperatures[i+1])
|
||||
|
||||
"""
|
||||
fig, (ax1) = P.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()
|
||||
"""
|
||||
|
||||
# Second step: Calcualte and Fit TNMH-model
|
||||
try:
|
||||
rate = float(str((re.findall(r'\d+\.\d+K', dataName))[0])[:-1])
|
||||
print("\tRate: %fK-min" % rate)
|
||||
except IndexError:
|
||||
rate = int(input("\tPlease enter the heating rate: "))
|
||||
print("\tRate: %fK-min" %rate)
|
||||
|
||||
|
||||
def tau_k(tau_0, deltaE, T_k, xx, T_f_km1):
|
||||
return tau_0 * np.exp(xx*deltaE/(R*T_k) + (1-xx)*deltaE/(R*T_f_km1))
|
||||
|
||||
|
||||
exponents = []
|
||||
taus = []
|
||||
temps = [0]
|
||||
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
|
||||
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 = [coordsLq[1][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)
|
||||
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]
|
||||
|
||||
print(T_f_ns[:10])
|
||||
return T_f_ns
|
||||
|
||||
|
||||
####################################################
|
||||
# Use T_g = T_fictive for tau determination
|
||||
T_g = T_f_0[0]
|
||||
# or plug in T_half, T_onset, ... from separate evaluation
|
||||
|
||||
# Start values
|
||||
p0 = [10, 0.5, 0.4, 225085]
|
||||
|
||||
|
||||
def fitTNMH(p, x):
|
||||
dTNMHdT = []
|
||||
modelTemp = np.append(x[::-1], x[1:])
|
||||
tau_0, xx, beta, deltaE = p
|
||||
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:]
|
||||
return np.array(res)
|
||||
|
||||
|
||||
# Use only every 20th data point to reduce fitting time
|
||||
temperaturesInterpol = np.linspace(temperaturesCut[0], temperaturesCut[len(temperaturesCut)-1], 200)
|
||||
dTfdTInterpol = np.interp(temperaturesInterpol, temperaturesCut, dTfdT)
|
||||
|
||||
ts = time.time()
|
||||
model = odr.Model(fitTNMH)
|
||||
whichFit = [1, 1, 1, 1]
|
||||
odrdata = odr.Data(temperaturesInterpol, dTfdTInterpol)
|
||||
myodr = odr.ODR(odrdata, model, beta0=p0, ifixb=whichFit, maxit=3000)
|
||||
result = (myodr.run())
|
||||
print(result.stopreason)
|
||||
p0 = result.beta
|
||||
print('TNMH heating Fit:')
|
||||
print(p0)
|
||||
# print('errors:')
|
||||
# print(result.sd_beta)
|
||||
ts = time.time()-ts
|
||||
print('Time for fit:')
|
||||
print(ts)
|
||||
tau_0, xx, beta, deltaE = p0
|
||||
|
||||
fitResult = fitTNMH(p0, temperaturesInterpol)
|
||||
|
||||
"""
|
||||
fig, (ax1) = P.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, '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))
|
||||
|
||||
save = 0 # input("Save parameter to file? (yes = 1, no = 0): ")
|
||||
if save == '1':
|
||||
np.savetxt(dataOutName, dataOut, fmt='%1.5e', header="tau_TNMH="+str(tau_0) + "\tbeta="+str(beta) + "\tx="+str(xx) + "\tdeltaE="+str(deltaE) + "\nTemp\t\tdTfdT_data\t\tdTfdT_fit")
|
||||
# Try to determine the heating rate from the data file name (Searching for digits followed by ".")
|
||||
paramName = ((str(dataName)).split("_"))[0] + "_DSC_tau.param"
|
||||
print(paramName)
|
||||
with open(paramName, 'a') as f:
|
||||
f.write('\n#Heating Rate \t tau_TNMH \t beta \t\t\t x \t deltaE \t Tg fictive')
|
||||
f.write('\n')
|
||||
f.write('%1.4e \t %1.4e \t %1.4e \t %1.4e \t %1.4e \t %1.4e' % (rate, tau_0, beta, xx, deltaE, T_g))
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user