From ec4094cd9238f695c1063efd6ab9a2b89005ee39 Mon Sep 17 00:00:00 2001
From: robrobo <robin.horstmann@gmx.net>
Date: Thu, 17 Jul 2025 12:34:43 +0200
Subject: [PATCH] added quick estimation for peak tau and height of non
 Gaussian

---
 src/mdevaluate/utils.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/src/mdevaluate/utils.py b/src/mdevaluate/utils.py
index 954590f..9bdbf2e 100644
--- a/src/mdevaluate/utils.py
+++ b/src/mdevaluate/utils.py
@@ -357,6 +357,38 @@ def quick1etau(t: ArrayLike, C: ArrayLike, n: int = 7) -> float:
     return tau_est
 
 
+def quicknongaussfit(t, C, width=4):
+    """
+    Estimates the time and height of the peak in the non-Gaussian function.
+    C is C(t) the correlation function
+    """
+    # TODO this is a very experimental interpolation, can fail
+    def ppoly(x,a,b,c,d,e,A,mu,sig):
+        return A*np.exp(-(x - mu)**2 / (2 * sig**2))+a+(b*x+e)*1/(1+np.exp(c*(x-d)))
+    # first rough estimate, the closest time. This is returned if the interpolation fails!
+    tau_est = t[np.argmax(C)]
+    nG_max = np.amax(C)
+    try:
+        with np.errstate(invalid='ignore'):
+            corr = C[t > 0]
+            time = np.log10(t[t > 0])
+            tau = time[np.argmax(corr)]
+            mask = (time>tau-width/2) & (time<tau+width/2)
+            time = time[mask] ; corr = corr[mask]
+            guess = [0.001,-0.001,5,tau-0.5,1.0,nG_max, tau, 1.5]
+            popt = curve_fit(ppoly, time, corr, p0=guess, maxfev=10000)[0]
+            # TODO instead use some root or solve
+            xspace = np.linspace(*time[[0,-1]], 10000)
+            y = ppoly(xspace, *popt)
+            tau_est = xspace[np.argmax(y)]
+            nG_max = np.amax(y)
+    except:
+        pass
+    if np.isnan(tau_est):
+        tau_est = np.inf
+    return tau_est, nG_max
+
+
 def susceptibility(
     time: NDArray, correlation: NDArray, **kwargs
 ) -> tuple[NDArray, NDArray]: