Closedown

script/motor-slide.py

@@ -0,0 +1,56 @@
+#Script imported from: PO2DV-NCS-LS_mot.xml
+
+#Pre-actions
+caput('PO2DV-NCS-LS:MOTOR.TWF', '0')
+sleep(0.5)
+caput('PO2DV-NCS-LS:MOTOR.RDBD', '0.1')
+
+#TODO: Set the diplay names of positioners and detectors
+scan = ManualScan(['VAL'], ['time', 'RVAL', 'TWF', 'RBV', 'Busy'] , [40.0], [44.0], [22])
+scan.start()
+
+#Creating channels: dimension 1
+#RegionPositioner VAL
+VAL = Channel('PO2DV-NCS-LS:MOTOR.VAL', type = 'd')
+VALReadback = Channel('PO2DV-NCS-LS:MOTOR.RBV', type = 'd')
+#Timestamp time
+#ScalarDetector RVAL
+RVAL = Channel('PO2DV-NCS-LS:MOTOR.RVAL', type = 'd')
+#ScalarDetector TWF
+TWF = Channel('PO2DV-NCS-LS:MOTOR.TWF', type = 'd')
+#ScalarDetector RBV
+RBV = Channel('PO2DV-NCS-LS:MOTOR.RBV', type = 'd')
+#ScalarDetector Busy
+Busy = Channel('PO2DV-NCS-LS:MOTOR.DMOV', type = 'd')
+
+#Dimension 1
+#RegionPositioner VAL
+for setpoint1 in frange(40.0, 42.0, 0.2, True) + frange(41.8, 40.0, 0.2, True):
+    VAL.put(setpoint1, timeout=None) # TODO: Set appropriate timeout
+    readback1 = VALReadback.get()
+    if abs(readback1 - setpoint1) > 0.1 : # TODO: Check accuracy
+        raise Exception('Actor VAL could not be set to the value ' + str(setpoint1))
+    #Detector time
+    detector1 = float(java.lang.System.currentTimeMillis())
+    #Detector RVAL
+    detector2 = RVAL.get()
+    #Detector TWF
+    detector3 = TWF.get()
+    #Detector RBV
+    detector4 = RBV.get()
+    #Detector Busy
+    detector5 = Busy.get()
+    scan.append ([setpoint1], [readback1], [detector1, detector2, detector3, detector4, detector5])
+
+#Closing channels
+VAL.close()
+VALReadback.close()
+RVAL.close()
+TWF.close()
+RBV.close()
+Busy.close()
+
+scan.end()
+
+#Post-actions
+caput('PO2DV-NCS-LS:MOTOR.RDBD', '1')

script/power-supply.py

@@ -0,0 +1,78 @@
+#Script imported from: PO2DV-NCS-HW_ps.xml
+from mathutils import estimate_peak_indexes, fit_gaussians, create_fit_point_list
+
+#Pre-actions
+caput('PO2DV-NCS-VHQ1:Set-RampA', '10')
+sleep(0.1)
+
+#TODO: Set the diplay names of positioners and detectors
+#ManualScan(writables, readables, start = None, end = None, steps = None, relative = False)
+scan = ManualScan(['time'], ['SetVA', 'ActualVA', 'ActualIA'] , [0.0], [20.0], [10]) #????????????????????? what does this do? what is writeables? what is readables?
+scan.start()
+
+#Creating channels: dimension 1
+#LinearPositioner SetVA
+SetVA = Channel('PO2DV-NCS-VHQ1:Set-VA', type = 'd')
+#Timestamp time
+#ScalarDetector ActualVA
+ActualVA = Channel('PO2DV-NCS-VHQ1:Actual-VA', type = 'd')
+#ScalarDetector ActualIA
+ActualIA = Channel('PO2DV-NCS-VHQ1:Actual-IA', type = 'd')
+
+#set voltage to 0
+print 'Ramping down power supply to 0V'
+SetVA.put(0.0, timeout=None)
+
+#wait up to 2 minutes for voltage to be ~0
+for setpoint1 in frange(0.0, 120.0, 1.0, True):
+    detector2 = ActualVA.get()
+    if detector2 <= 1.0:
+        break
+    sleep(0.5)
+
+#Dimension 1
+#LinearPositioner SetVA
+print 'Ramping up power supply'
+for setpoint1 in frange(0.0, 20.0, 10.0, True): #?????????????????? What is the relationship between this FOR loop and the ManualScan?
+    if setpoint1 > 50.0 or setpoint1 < 0.0:
+        break
+    SetVA.put(setpoint1, timeout=None) # TODO: Set appropriate timeout
+    readback1 = SetVA.get()
+    if abs(readback1 - setpoint1) > 0.5 : # TODO: Check accuracy
+        raise Exception('Actor SetVA could not be set to the value ' + str(setpoint1))
+    #scan quickly the output during some seconds
+    for setpoint2 in range(0, 20):
+        #Detector time
+        detector1 = float(java.lang.System.currentTimeMillis())
+        #Detector ActualVA
+        detector2 = ActualVA.get()
+        detector3 = ActualIA.get()
+        #scan.append ([setpoint1], [readback1], [detector1, detector2])
+
+        #append(setpoints, positions, values)
+        scan.append ([detector1], [detector1], [readback1, detector2, detector3]) #?????????????????????? what is setpoint? is position = X-axis?
+        sleep( 0.1 ) # Settling time
+
+#reset output to 0V
+SetVA.put(0.0, timeout=None)
+#Closing channels
+SetVA.close()
+ActualVA.close()
+ActualIA.close()
+
+readable =    scan.readables[0]
+positions =   scan.writables
+
+threshold = (min(readable) + max(readable))/2
+min_peak_distance = 5.0
+
+peaks = estimate_peak_indexes(readable, positions, threshold, min_peak_distance)
+print "Peak indexes: " + str(peaks)
+print "Peak x: " + str(map(lambda x:positions[x], peaks))
+print "Peak y: " + str(map(lambda x:readable[x], peaks))
+
+
+scan.end()
+
+
+#??????????????????????????? Device modeling: how can I modelise a Power Supply, for example?

script/test/data.py

@@ -0,0 +1,57 @@
+data = [1,2,3,4,5]
+path="group/data"
+
+
+save_dataset(path, data)
+read =load_data(path)
+
+print read.tolist()
+
+#plot(read)
+
+
+data = [ [1,2,3,4,5],  [2,3,4,5,6],  [3,4,5,6,7]]
+path="group/data2"
+
+
+save_dataset(path, data)
+set_attribute(path, "AttrString", "Value")
+set_attribute(path, "AttrInteger", 1)
+set_attribute(path, "AttrDouble", 2.0)
+set_attribute(path, "AttrBoolean", True)
+
+read =load_data(path)
+
+print read.tolist()
+
+plot(read)
+
+path = "group/data3"
+create_dataset(path, 'i')
+for i in range(10):
+    save_data_item(path,i)    
+
+
+path = "group/data4"
+create_dataset(path, 'd', False, (0,0))
+for row in data:
+    save_data_item(path, row)    
+
+path = "group/data5"
+names = ["a", "b", "c", "d"]
+types = ["d", "d", "d", "[d"]
+lenghts = [0,0,0,5]
+dims = [0,0,0,0]
+
+
+data = [    [1,2,3,[0,1,2,3,4]],  
+            [2,3,4,[3,4,5,6,7]],  
+            [3,4,5,[6,7,8,9,4]]    ]
+
+create_table(path, names, types, lenghts, dims)
+for row in data:
+    save_table_item(path, row)   
+flush_data() 
+
+read =load_data(path)
+print read

script/test/import.py

@@ -0,0 +1,42 @@
+#Script imported from: test1.xml
+
+#Variables
+var1 = 0.0
+
+#TODO: Set the diplay names of positioners and detectors
+scan = ManualScan(['id278043'], ['id348623', 'id367393'] , [0.0], [31.0], [31])
+scan.start()
+
+#Creating channels: dimension 1
+#LinearPositioner id278043
+id278043 = Channel('TESTIOC:TESTCALCOUT:Input', type = 'd')
+#ScalarDetector id348623
+id348623 = Channel('TESTIOC:TESTCALCOUT:Output', type = 'd')
+#ScalarDetector id367393
+id367393 = Channel('TESTIOC:TESTSINUS:SinCalc', type = 'd')
+
+#Dimension 1
+#LinearPositioner id278043
+for setpoint1 in frange(0.0, 31.0, 1.0, True):
+    if setpoint1 > 31.0 or  setpoint1 < 0.0:
+        break
+    id278043.put(setpoint1, timeout=None) # TODO: Set appropriate timeout
+    readback1 = id278043.get()
+    if abs(readback1 - setpoint1) > 0.5 : # TODO: Check accuracy
+        raise Exception('Actor id278043 could not be set to the value ' + str(setpoint1))
+    #Dimension Actions
+    #Script action
+    #TODO: Move, if needed, this import to the file header: import time
+    time.sleep(0.1)
+    #Detector id348623
+    detector1 = id348623.get()
+    #Detector id367393
+    detector2 = id367393.get()
+    scan.append ([setpoint1], [readback1], [detector1, detector2])
+
+#Closing channels
+id278043.close()
+id348623.close()
+id367393.close()
+
+scan.end()

script/test/import2.py

@@ -0,0 +1,55 @@
+#Script imported from: test2.xml
+
+#Variables
+var1 = 0.0
+
+#TODO: Set the diplay names of positioners and detectors
+scan = ManualScan(['id643271', 'id278043'], ['id348623', 'id367393'] , [0.0, 0.0], [5.0, 31.0], [5, 31])
+scan.start()
+
+#Creating channels: dimension 1
+#LinearPositioner id643271
+id643271 = Channel('TESTIOC:TESTCALC:MyCalc', type = 'd')
+#Creating channels: dimension 2
+#LinearPositioner id278043
+id278043 = Channel('TESTIOC:TESTCALCOUT:Input', type = 'd')
+#ScalarDetector id348623
+id348623 = Channel('TESTIOC:TESTCALCOUT:Output', type = 'd')
+#ScalarDetector id367393
+id367393 = Channel('TESTIOC:TESTSINUS:SinCalc', type = 'd')
+
+#Dimension 1
+#LinearPositioner id643271
+for setpoint1 in frange(0.0, 5.0, 1.0, True):
+    if setpoint1 > 5.0 or  setpoint1 < 0.0:
+        break
+    id643271.put(setpoint1, timeout=None) # TODO: Set appropriate timeout
+    readback1 = id643271.get()
+    if abs(readback1 - setpoint1) > 0.5 : # TODO: Check accuracy
+        raise Exception('Actor id643271 could not be set to the value ' + str(setpoint1))
+    #Dimension 2
+    #LinearPositioner id278043
+    for setpoint2 in frange(0.0, 31.0, 1.0, True):
+        if setpoint2 > 31.0 or  setpoint2 < 0.0:
+            break
+        id278043.put(setpoint2, timeout=None) # TODO: Set appropriate timeout
+        readback2 = id278043.get()
+        if abs(readback2 - setpoint2) > 0.5 : # TODO: Check accuracy
+            raise Exception('Actor id278043 could not be set to the value ' + str(setpoint2))
+        #Dimension Actions
+        #Script action
+        #TODO: Move, if needed, this import to the file header: import time
+        time.sleep(0.01)
+        #Detector id348623
+        detector1 = id348623.get()
+        #Detector id367393
+        detector2 = id367393.get()
+        scan.append ([setpoint1, setpoint2], [readback1, readback2], [detector1, detector2])
+
+#Closing channels
+id278043.close()
+id348623.close()
+id367393.close()
+id643271.close()
+
+scan.end()

script/test/import5.py

@@ -0,0 +1,31 @@
+#Script imported from: test5.xml
+
+set_flush_records(False)
+#set_preference("scanPlotDisabled",True)
+set_preference("scanTableDisabled",True)
+
+#TODO: Set the diplay names of positioners and detectors
+scan = ManualScan(['id236750', 'id226053'], ['id246209'] , [0.0, 0.0], [1000.0, 1000.0], [1000, 1000])
+scan.start()
+
+#Creating channels: dimension 1
+#PseudoPositioner id236750
+#Timestamp id246209
+#Creating channels: dimension 2
+#PseudoPositioner id226053
+
+#Dimension 1
+#PseudoPositioner id236750
+for setpoint1 in range(0, 1000):
+    readback1 = setpoint1
+    #Detector id246209
+    detector1 = java.lang.System.currentTimeMillis()
+    #Dimension 2
+    #PseudoPositioner id226053
+    for setpoint2 in range(0, 1000):
+        readback2 = setpoint2
+        scan.append ([setpoint1, setpoint2], [readback1, readback2], [detector1])
+
+#Closing channels
+
+scan.end()

script/test/import7.py

@@ -0,0 +1,52 @@
+#Script imported from: test7.xml
+
+#Variables
+var1 = 0.0
+
+#TODO: Set the diplay names of positioners and detectors
+scan = ManualScan(['id278043'], ['id348623', 'id367393', 'id192931', 'id173831'] , [0.0], [31.0], [31])
+scan.start()
+
+#Creating channels: dimension 1
+#LinearPositioner id278043
+id278043 = Channel('TESTIOC:TESTCALCOUT:Input', type = 'd')
+#ScalarDetector id348623
+id348623 = Channel('TESTIOC:TESTCALCOUT:Output', type = 'd')
+#ScalarDetector id367393
+id367393 = Channel('TESTIOC:TESTSINUS:SinCalc', type = 'd')
+#Timestamp id192931
+
+#Dimension 1
+#LinearPositioner id278043
+for setpoint1 in frange(0.0, 31.0, 1.0, True):
+    if setpoint1 > 31.0 or  setpoint1 < 0.0:
+        break
+    id278043.put(setpoint1, timeout=None) # TODO: Set appropriate timeout
+    readback1 = id278043.get()
+    if abs(readback1 - setpoint1) > 0.5 : # TODO: Check accuracy
+        raise Exception('Actor id278043 could not be set to the value ' + str(setpoint1))
+    #Dimension Actions
+    #Script action
+    #Variable Mappings
+    x = Channel('TESTIOC:TESTSINUS:SinCalc', type = 'd')
+    #TODO: Move, if needed, this import to the file header: import time
+    time.sleep(0.1)
+    #print "==>" + str(x.getValue())
+    #Detector id348623
+    detector1 = id348623.get()
+    #Detector id367393
+    detector2 = id367393.get()
+    #Detector id192931
+    detector3 = float(java.lang.System.currentTimeMillis())
+    #Manipulation id173831
+    #Variable Mappings
+    v =  detector2
+    id173831 = v+100.0
+    scan.append ([setpoint1], [readback1], [detector1, detector2, detector3, id173831])
+
+#Closing channels
+id278043.close()
+id348623.close()
+id367393.close()
+
+scan.end()

script/test/parallel.py

@@ -0,0 +1,42 @@
+
+#Simple parallization
+def task1():
+    return out.read()
+
+def task2():
+    return inp.read()
+
+def task3():
+    time.sleep(0.1)
+    return sin.read()
+
+ret = parallelize(task1, task2, task3)
+print ret
+
+
+#Fork amd join
+ret = fork(task1, task2, task3)
+time.sleep(0.1)
+ret = join(ret)
+print ret
+
+
+#Functions with parameters
+def devRead(dev, msg):
+    print msg + " -> " + dev.getName()
+    return dev.read()
+
+ret = parallelize((devRead,(out,"1")), (devRead,(inp,"2")), (devRead,(sin,"3")))
+print ret
+
+
+#Exception in parallel task
+def taskExcept(msg):
+    raise Exception ("Error in parallel task " + msg)
+
+
+ret = parallelize((taskExcept,("1")), (taskExcept,(inp,"2")) )
+print ret
+
+
+

script/test/scan.py

@@ -0,0 +1,3 @@
+def calc(a):
+    return a*2
+

script/test/test0.py

@@ -0,0 +1,60 @@
+import sys
+import time
+
+#To add library folders from within the script
+#sys.path.append("./site-packages")
+
+
+import requests
+r = requests.get('https://api.github.com', auth=('user', 'pass')) 
+print r.status_code
+print r.headers['content-type']
+r.close()
+
+
+def calc2(a):
+    return a*2
+time.sleep(0.1)
+#import os
+#print os.environ
+#import calc
+
+
+time.sleep(0.1)
+for x in range(3):
+    print x
+    while(True):
+        print x*2
+        break
+time.sleep(2)
+
+
+data = [1,2,3,4,5]
+path="group/data"
+
+
+save_dataset(path, data)
+read =load_data(path)
+
+print read.tolist()
+
+#plot(read)
+
+
+data = [ [1,2,3,4,5],  [2,3,4,5,6],  [3,4,5,6,7]]
+path="group/data2"
+
+
+save_dataset(path, data)
+read =load_data(path)
+
+print read.tolist()
+
+plot(read)
+
+
+
+
+"""
+It lives!!!!
+"""

script/test/test1.py

@@ -0,0 +1,15 @@
+""" 
+Line Scan 
+"""
+#setPreference("enabledPlots",(sin,out))
+
+#def lscan(writables, readables, start, end, steps, latency=0.0, relative = False, context=None, before_read=None, after_read=None):
+a= lscan(inp, (sin,out,arr), 0, 40, 10, 0.1, context = "Scan1")
+
+path = get_current_data_group()
+set_attribute(path, "AttrString", "Value")
+set_attribute(path, "AttrInteger", 1)
+set_attribute(path, "AttrDouble", 2.0)
+set_attribute(path, "AttrBoolean", True)
+
+

script/test/test10.py

@@ -0,0 +1,44 @@
+"""
+Multi-peak search 
+""" 
+
+from mathutils import estimate_peak_indexes, fit_gaussians, create_fit_point_list
+
+start = 0
+end = 30
+step_size = 0.2
+
+result= lscan(sout,sinp,start,end,[step_size,],0.01)
+
+readable = result.getReadable(0)
+positions = result.getPositions(0)
+
+threshold = (min(readable) + max(readable))/2
+min_peak_distance = 5.0
+
+peaks = estimate_peak_indexes(readable, positions, threshold, min_peak_distance)
+print "Peak indexes: " + str(peaks)
+print "Peak x: " + str(map(lambda x:positions[x], peaks))
+print "Peak y: " + str(map(lambda x:readable[x], peaks))
+
+
+
+gaussians = fit_gaussians(readable, positions, peaks)
+
+
+plots = plot([readable],["sin"],[positions] )
+for i in range(len(peaks)):
+    peak = peaks[i]
+    (norm, mean, sigma) = gaussians[i]
+    print 'gaussians: ' + str(gaussians[i])
+    if abs(mean - positions[peak]) < min_peak_distance:
+        print "Peak -> " +  str(mean)
+        plots[0].addMarker(mean, None, "N="+str(round(norm,2)), None)
+    else:
+        print "Invalid gaussian fit: " +  str(mean)
+
+print "OK
+
+
+
+

script/test/test11.py

@@ -0,0 +1,11 @@
+""" 
+Parameters
+"""
+
+start = float(args[0])
+end = float(args[1])
+step =  int(args[2])
+
+a= lscan(inp, sin, start, end, step, 0.1)
+
+

script/test/test12.py

@@ -0,0 +1,8 @@
+""" 
+Relative Scan 
+"""
+
+a= lscan(inp, (sin,out,arr), -2, 2, 20, 0.1, True)
+
+
+

script/test/test13.py

@@ -0,0 +1,8 @@
+""" 
+Relative Scan 
+"""
+
+a= lscan(inp, (sin,out,arr), -2, 2, 20, 0.1, True)
+
+
+

script/test/test2.py

@@ -0,0 +1,16 @@
+""" 
+Line Scan with 2 writables and triggering
+"""
+
+index=0 
+
+def BeforeReadout():
+    global index
+    print "Frame = " + str(index)
+    index=index+1
+    #log("trigger " + index)
+    caput("TESTIOC:TESTBO:MyBO","On")
+    caput("TESTIOC:TESTBO:MyBO","Off")    
+
+
+a= lscan((motor,inp), (sin,out,arr), (0,0), (4,40), 20, 0.1, before_read=BeforeReadout)

script/test/test3.py

@@ -0,0 +1,21 @@
+"""
+Processing and plotting scan data
+""" 
+
+inp.write(0.0)
+scan1= lscan(inp, (sin,out,arr), 0, 40, 20, 0.1, False, "Scan 1")
+scan2= lscan(inp, (sin,out,arr), 0, 40, 20, 0.1, False, "Scan 2")
+
+
+from operator import add
+result = map(add, scan1.getReadable(0), scan2.getReadable(0))
+
+#Alternative:
+#result=[]
+#for i in range(len(scan1.records)):
+#	result.append(scan1.records[i].values[0]+scan2.records[i].values[0])
+
+
+plot(result)
+print result
+

script/test/test4.py

@@ -0,0 +1,12 @@
+"""
+Vector Scan
+"""
+
+vector = [  [1,1] , [1,2] , [1,3] , [1,4]   ,
+                     [1.5,2.5]              , 
+            [2,1] , [2,2] , [2,3] , [2,4]   ,
+                     [2.5,2.5]              ,
+            [3,1] , [3,2] , [3,3] , [3,4]     ]
+
+#def vscan(writables, readables, vector, latency=0.0, relative = False, context=None, before_read=None, after_read=None):
+a= vscan((dev,inp),(sin,out),vector,0.5)

script/test/test5.py

@@ -0,0 +1,5 @@
+"""
+Area Scan
+"""
+
+ascan((dev,out), (sin,out,arr), (0,10), (20,30), (20,20))

script/test/test6.py

@@ -0,0 +1,36 @@
+""" 
+Creating pseudo-devices
+"""
+import time
+
+sin_val=None        
+        
+class Sensor(ch.psi.pshell.dev.Readable):
+    def read(self):
+        global sin_val
+        return sin_val + time.clock()
+    def getName(self):
+        return "Sensor"
+
+class Positioner(ch.psi.pshell.dev.Writable):
+    def write(self,pos):
+        print pos
+    def getName(self):
+        return "Positioner"
+
+class Listener (ch.psi.pshell.dev.DeviceListener):
+    def onStateChanged(self, device, state, former):
+        pass
+    def onValueChanged(self, device, value, former):
+        global sin_val
+        sin_val=value
+
+sensor=Sensor()
+positioner=Positioner()
+listener = Listener()
+
+sin.addListener(listener)
+try:
+    a= lscan((inp,positioner),(sin,sensor),(0,0),(40,10),20,0.1)
+finally:
+    sin.removeListener(listener)

script/test/test7.py

@@ -0,0 +1,114 @@
+"""
+Function fitting and peak search directly with org.apache.commons.math3
+""" 
+
+import org.apache.commons.math3.util.FastMath as FastMath
+
+import org.apache.commons.math3.stat.regression.SimpleRegression as SimpleRegression
+
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction as PolynomialFunction
+import org.apache.commons.math3.analysis.function.Gaussian as Gaussian
+import org.apache.commons.math3.analysis.function.HarmonicOscillator as HarmonicOscillator 
+import org.apache.commons.math3.analysis.solvers.LaguerreSolver as LaguerreSolver
+
+import org.apache.commons.math3.fitting.GaussianCurveFitter as GaussianCurveFitter
+import org.apache.commons.math3.fitting.PolynomialCurveFitter as PolynomialCurveFitter
+import org.apache.commons.math3.fitting.HarmonicCurveFitter as HarmonicCurveFitter
+import org.apache.commons.math3.fitting.WeightedObservedPoint as WeightedObservedPoint
+
+import org.apache.commons.math3.analysis.differentiation.FiniteDifferencesDifferentiator as FiniteDifferencesDifferentiator
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure as DerivativeStructure
+
+ 
+start = 0
+end = 2
+step_size = 0.1
+
+result= lscan(inp,(sin,out),start,end,[step_size,],0.1)
+
+readable = result.getReadable(0)
+positions = result.getPositions(0)
+
+simple_regression= SimpleRegression()
+
+values = []
+for record in result.records:
+    simple_regression.addData(record.positions[0], record.values[0])
+    values.append(WeightedObservedPoint(1, record.positions[0], record.values[0]))
+
+order = 6
+#polynomial_fitter =  PolynomialCurveFitter.create(0).withStartPoint([ 5.0, 5.0, 5.0])
+polynomial_fitter = PolynomialCurveFitter.create(order)
+gaussian_fitter = GaussianCurveFitter.create();
+harmonic_fitter = HarmonicCurveFitter.create();
+
+
+best = polynomial_fitter.fit(values)
+fitted_polynomial_function = PolynomialFunction(best)
+
+
+(normalization, mean, sigma) = gaussian_fitter.fit(values)
+print (normalization, mean, sigma) 
+fitted_gaussian_function = Gaussian(normalization, mean, sigma)
+print "Mean: " + str(mean)
+
+(amplitude, angular_frequency, phase) = harmonic_fitter.fit(values)
+print (amplitude, angular_frequency, phase) 
+fitted_harmonic_function = HarmonicOscillator(amplitude, angular_frequency, phase)
+
+
+#differentiator = FiniteDifferencesDifferentiator(order+2, 0.25) #points, step size
+#derivative = differentiator.differentiate(fitted_polynomial)
+derivative = fitted_polynomial_function.derivative() 
+print fitted_polynomial_function.coefficients 
+print derivative.coefficients 
+
+
+regression = []
+fit_polinomial = []
+fit_gaussian = []
+fit_harmonic = []
+fit_polinomial_derivative = []
+resolution = step_size/100
+for x in frange(start,end,resolution, True):
+    regression.append(simple_regression.predict(x))
+    fit_polinomial.append(fitted_polynomial_function.value(x))
+    fit_gaussian.append(fitted_gaussian_function.value(x))
+    fit_harmonic.append(fitted_harmonic_function.value(x))
+    fit_polinomial_derivative.append(derivative.value(x))
+x = frange(start, end+resolution, resolution)
+
+
+
+
+def calculate_peaks(function, start_value = sys.float_info.min, end_value = sys.float_info.max, positive=True):
+    derivative = function.derivative() 
+    derivative2 = derivative.derivative() 
+    ret = []
+    solver = LaguerreSolver()
+    for complex in solver.solveAllComplex(derivative.coefficients, start_value):
+        r = complex.real
+        if start_value < r < end_value:
+            if (positive and (derivative2.value(r) < 0)) or ( (not positive) and (derivative2.value(r) > 0)):         
+                ret.append(r)
+    return ret
+
+peaks = calculate_peaks(fitted_polynomial_function)
+
+
+plots = plot([readable, regression, fit_polinomial, fit_gaussian, fit_harmonic, fit_polinomial_derivative] , 
+    ["sin", "regression", "fit polinomial", "fit gaussian", "fit harmonic ", "fit polinomial derivative"], 
+    xdata = [positions,x,x,x,x,x] )
+
+
+for p in peaks:    
+    print "Max: " + str(p)
+    plots[0].addMarker(p, None, "Max", None)
+
+#plots[0].addMarker(mean, None, "Gaussian Mean", None)
+
+
+
+print result
+
+

script/test/test8.py

@@ -0,0 +1,36 @@
+"""
+Multi-peak search 
+""" 
+
+from mathutils import estimate_peak_indexes, fit_gaussians
+
+start = 0
+end = 5
+step_size = 0.1
+
+result= lscan(inp,sin,start,end,[step_size,],0.1)
+
+readable = result.getReadable(0)
+positions = result.getPositions(0)
+
+threshold = (min(readable) + max(readable))/2
+min_peak_distance = 0.5
+
+peaks = estimate_peak_indexes(readable, positions, threshold, min_peak_distance)
+print peaks
+
+gaussians = fit_gaussians(readable, positions, peaks)
+
+plots = plot([readable],["sin"],[positions] )
+for i in range(len(peaks)):
+    peak = peaks[i]
+    (norm, mean, sigma) = gaussians[i]
+    if abs(mean - positions[peak]) < min_peak_distance:
+        print "Peak -> " +  str(mean)
+        plots[0].addMarker(mean, None, "N="+str(round(norm,2)), None)
+    else:
+        print "Invalid gaussian fit: " +  str(mean)
+
+
+
+

script/test/test9.py

@@ -0,0 +1,58 @@
+"""
+Function fitting and peak search with mathutils facade
+""" 
+from mathutils import fit_polynomial,fit_gaussian, fit_harmonic, calculate_peaks
+from mathutils import PolynomialFunction, Gaussian, HarmonicOscillator
+import math
+
+
+start = 0
+end = 10
+step_size = 0.1
+
+result= lscan(sout,sinp,start,end,[step_size,],0.01)
+
+readable = result.getReadable(0)
+positions = result.getPositions(0)
+
+def get_function_data(function, start, end, resolution):
+    ret = []
+    for x in frange(start, end, resolution, True):
+        fit_polinomial.append(function.value(x))
+
+
+pars_polynomial = (a0, a1, a2, a3, a4, a5, a6) = fit_polynomial(readable, positions, 6)
+fitted_polynomial_function = PolynomialFunction(pars_polynomial)
+print pars_polynomial 
+
+(normalization, mean, sigma) = fit_gaussian(readable, positions, True)
+fitted_gaussian_function = Gaussian(normalization, mean, sigma)
+print (normalization, mean, sigma)
+
+(amplitude, angular_frequency, phase) = fit_harmonic(readable, positions)
+fitted_harmonic_function = HarmonicOscillator(amplitude, angular_frequency, phase)
+print (amplitude, angular_frequency, phase) 
+
+
+resolution = step_size/100
+fit_polinomial = []
+fit_gaussian = []
+fit_harmonic = []
+for x in frange(start,end,resolution, True):
+    fit_polinomial.append(fitted_polynomial_function.value(x))
+    fit_gaussian.append(fitted_gaussian_function.value(x))
+    fit_harmonic.append(fitted_harmonic_function.value(x))
+x = frange(start, end+resolution, resolution)
+
+
+peaks = calculate_peaks(fitted_polynomial_function)
+
+plots = plot([readable, fit_polinomial, fit_gaussian, fit_harmonic] , 
+    ["data", "polinomial", "gaussian", "harmonic"], xdata = [positions,x,x,x] )
+
+for p in peaks:    
+    print "Max: " + str(p)
+    plots[0].addMarker(p, None, "Max=" + str(round(p,2)), None)
+import java.awt.Color
+plots[0].addMarker(mean, None, "Mean=" + str(round(mean,2)), java.awt.Color.LIGHT_GRAY)
+

script/test/testall.py

@@ -0,0 +1,6 @@
+
+scripts = ["test11", "test1", "test2", "test3", "test4", "test5", "test6",  "test7", "test8", "test9", "test10", "import", "import2", "import7", "data", "test0"]
+
+for script in scripts:
+    set_status("Running: " + script)
+    run(script)