Startup

script/CPython/gfitoff.py

@@ -0,0 +1,30 @@
+import numpy as np
+import scipy.optimize
+
+
+def gfitoff(x, y, off=None, amp=None, com=None, sigma=None):
+    if off is None:
+        off = y.min()  # good enough starting point for offset
+
+    if com is None:
+        com = x[y.argmax()]
+
+    if amp is None:
+        amp = y.max() - off
+
+    # For normalised gauss curve sigma=1/(amp*sqrt(2*pi))
+    if sigma is None:
+        surface = np.trapz((y-off), x=x)
+        sigma = surface / (amp * np.sqrt(2 * np.pi))
+    try:
+        popt, pcov = scipy.optimize.curve_fit(gauss_fn, x, y, p0=[off, amp, com, sigma], method='lm')
+        popt[3] = abs(popt[3])  # sigma should be returned as positive
+    except Exception as e:
+        print("Gauss fitting not successful.\n" + str(e))
+        popt = [off, amp, com, abs(sigma)]
+
+    return popt
+
+
+def gauss_fn(x, a, b, c, d):
+    return a + b * np.exp(-(np.power((x - c), 2) / (2 * np.power(d, 2))))

script/CPython/hfitoff.py

@@ -0,0 +1,23 @@
+import numpy as np
+from scipy.optimize import leastsq
+
+
+def hfitoff(data, xdeg):
+    """
+    Harmonic fit with offset
+    """
+    start, end = min(xdeg), max(xdeg)
+    guess_amplitude = 2 ** 0.5 * np.std(data)
+    guess_phase = 0
+    guess_offset = np.mean(data)
+    xrad = xdeg / 180 * np.pi
+    optimize_func = lambda x: x[0] * np.sin(xrad + x[1]) + x[2] - data
+    fit_amplitude, fit_phase, fit_offset = leastsq(optimize_func, [guess_amplitude, guess_phase, guess_offset])[0]
+    fit_xdeg = np.linspace(start, end, 100)
+    fit_data = fit_amplitude * np.sin(fit_xdeg / 180 * np.pi + fit_phase) + fit_offset
+    ph_crest = 90 - fit_phase / np.pi * 180
+    xdeg_max = xdeg[np.argmax(data)]
+    if start <= ph_crest <= end:
+        return (fit_amplitude, fit_phase, fit_offset, True, ph_crest, fit_xdeg, fit_data)
+    else:
+        return (fit_amplitude, fit_phase, fit_offset, False, xdeg_max, fit_xdeg, fit_data)

script/CPython/wrapper.py

@@ -0,0 +1,10 @@
+from jeputils import *
+
+def hfitoff(data, xdeg):
+    ret =  call_jep("CPython/hfitoff", "hfitoff", [to_npa(data),to_npa(xdeg)])
+    print ret
+    return (ret[0], ret[1], ret[2],ret[3], ret[4],  ret[5].data,ret[6].data)
+
+def gfitoff(x, y, off=None, amp=None, com=None, sigma=None):
+    ret = call_jep("CPython/gfitoff", "gfitoff", [to_npa(x), to_npa(y), off, amp, com, sigma])
+    return ret if ret is None or is_list(ret) else ret.data