tg evaluation in dsc

src/nmreval/data/dsc.py

@@ -4,6 +4,7 @@ from typing import Optional
 
 import numpy as np
 from scipy.optimize import fsolve
+from scipy.signal import savgol_filter
 
 try:
     from scipy.integrate import cumulative_trapezoid
@@ -23,6 +24,11 @@ class DSC(Points):
 
         x, unique = np.unique(x, return_index=True)
 
+        if 'y_err' in kwargs:
+            _yerr = kwargs['y_err']
+            _yerr = np.asarray(_yerr).reshape(np.asarray(x).shape)
+            kwargs['y_err'] = _yerr[unique]
+
         super().__init__(x, y[unique], **kwargs)
 
     def get_fictive_cp(self, glass: tuple[float, float], liquid: tuple[float, float]) -> ('DSC', float):
@@ -44,17 +50,18 @@ class DSC(Points):
         real_area = cumulative_trapezoid(region.y, region.x, initial=0)
         real_area -= real_area[-1]
 
-        t = regress2.intercept-regress.intercept
+        t = regress2.intercept - regress.intercept
-        m = regress2.slope-regress.slope
+        m = regress2.slope - regress.slope
-        c0 = 0.5*m*region.x.max()**2 + t*region.x.max()
+        c0 = 0.5 * m * region.x.max() ** 2 + t * region.x.max()
 
         def equiv(_x, _i):
-            return (0.5*m * _x**2 + t*_x - c0) - real_area[_i]
+            return (0.5 * m * _x ** 2 + t * _x - c0) - real_area[_i]
 
         def equiv_prime(_x, _i):
             return m * _x + t
 
-        fictive_temperature = np.array([fsolve(equiv, region.x[i], fprime=equiv_prime, args=(i,))[0] for i in range(len(region.x))])
+        fictive_temperature = np.array(
+            [fsolve(equiv, region.x[i], fprime=equiv_prime, args=(i,))[0] for i in range(len(region.x))])
 
         t_g_fictive = fictive_temperature[:20].mean()
         region.y = np.gradient(fictive_temperature, region.x)
@@ -62,7 +69,8 @@ class DSC(Points):
         return region, t_g_fictive
 
     def calculate_tnmh(self, p0: list, glass: tuple[float, float], liquid: tuple[float, float],
-                       tg: float = None, num_points: int = 200, return_fictive: bool = True) -> ('FitResult', Optional[float], Optional[DSC]):
+                       tg: float = None, num_points: int = 200, return_fictive: bool = True) \
+            -> ('FitResult', Optional[float], Optional[DSC]):
 
         dtf_dt, fictive_tg = self.get_fictive_cp(glass, liquid)
         if tg is None:
@@ -75,7 +83,7 @@ class DSC(Points):
         fitter = FitRoutine()
         fitter.set_model(TNMH)
         data = fitter.add_data(temp_equidist, dtf_dt_equidist)
-        data.set_parameter(p0 + [tg, self.value], var=[True]*4 + [False]*2, default_bounds=True)
+        data.set_parameter(p0 + [tg, self.value], var=[True] * 4 + [False] * 2, default_bounds=True)
 
         res = fitter.run()[0]
 
@@ -83,3 +91,41 @@ class DSC(Points):
             return res, tg, dtf_dt
         else:
             return res
+
+    def glass_transition(self, glass, liquid):
+        low_idx = tuple(np.argmin(np.abs(self.x - g)) for g in glass)
+        high_idx = tuple(np.argmin(np.abs(self.x - l)) for l in liquid)
+
+        x = self.x[low_idx[0]:high_idx[1]]
+        y = self.y[low_idx[0]:high_idx[1]]
+
+        yy = savgol_filter(y, window_length=min(len(x) // 20, 50), polyorder=1, deriv=1) / np.mean(np.diff(x))
+
+        high_idx = (high_idx[0] - low_idx[0], high_idx[1] - low_idx[0])
+        low_idx = (0, low_idx[1] - low_idx[0])
+
+        inflection = np.argmax(yy)
+
+        p1 = linregress(x[low_idx[0]:low_idx[1]], y[low_idx[0]:low_idx[1]])
+        glass_baseline = p1.slope * x + p1.intercept
+
+        p2 = linregress(x[high_idx[0]:high_idx[1]], y[high_idx[0]:high_idx[1]])
+        liquid_baseline = p2.slope * x + p2.intercept
+
+        tangent_line = yy[inflection] * (x - x[inflection]) + y[inflection]
+
+        onset = np.argmin(np.abs(tangent_line - glass_baseline))
+        end = np.argmin(np.abs(tangent_line - liquid_baseline))
+
+        midpoint = np.argmin(np.abs(y - 0.5 * (liquid_baseline[end] - glass_baseline[onset])))
+        cut_tangent = np.where((tangent_line > y.min() - 1) & (tangent_line < y.max() + 1))
+
+        return {
+            'onset': (x[onset], glass_baseline[onset]),
+            'mid': (x[midpoint], y[midpoint]),
+            'end': (x[end], liquid_baseline[end]),
+            'inflection': (x[inflection], y[inflection]),
+            'glass_baseline': np.c_[x, glass_baseline],
+            'liquid_baseline': np.c_[x, liquid_baseline],
+            'tangent_line': np.c_[x[cut_tangent], tangent_line[cut_tangent]],
+        }