From ef32e7aa7a441b3e4790d0d6b7111f9ac866fe8c Mon Sep 17 00:00:00 2001
From: Sebastien Bourdeauducq <sb@m-labs.hk>
Date: Wed, 21 Jan 2015 14:35:37 +0800
Subject: [PATCH] examples/flopping_f_simulation: fitting

---
 examples/flopping_f_simulation.py | 27 +++++++++++++++++++--------
 examples/pdb.pyon                 |  2 +-
 2 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/examples/flopping_f_simulation.py b/examples/flopping_f_simulation.py
index 6f1de1567..e2b13a349 100644
--- a/examples/flopping_f_simulation.py
+++ b/examples/flopping_f_simulation.py
@@ -1,10 +1,18 @@
 from math import sqrt, cos, pi
 import time
+import random
+
+import numpy as np
+from scipy.optimize import curve_fit
 
 from artiq import *
 
 
-def model(x, F0=1500, A=80, B=40, t=0.02, tpi=0.03):
+def model(x, F0):
+    t = 0.02
+    tpi = 0.03
+    A = 80
+    B = 40
     return A+(B-A)/2/(4*tpi**2*(x-F0)**2+1)*(1-cos(pi*t/tpi*sqrt(4*tpi**2*(x-F0)**2+1)))
 
 
@@ -16,6 +24,9 @@ class FloppingF(AutoDB):
         min_freq = Argument(1000)
         max_freq = Argument(2000)
 
+        F0 = Argument(1500)
+        noise_amplitude = Argument(0.1)
+
         frequency = Result()
         brightness = Result()
 
@@ -30,16 +41,16 @@ class FloppingF(AutoDB):
     def run(self):
         for i in range(self.npoints):
             frequency = (self.max_freq-self.min_freq)*i/(self.npoints - 1) + self.min_freq
-            brightness = model(frequency)
+            brightness = model(frequency, self.F0) + self.noise_amplitude*random.random()
             self.frequency.append(frequency)
             self.brightness.append(brightness)
             time.sleep(0.1)
         self.analyze()
 
     def analyze(self):
-        min_f = self.frequency.read[0]
-        min_b = self.brightness.read[0]
-        for f, b in zip(self.frequency.read, self.brightness.read):
-            if b < min_b:
-                min_f, min_b = f, b
-        self.flopping_freq = min_f
+        popt, pcov = curve_fit(lambda xdata, F0: [model(x, F0) for x in xdata],
+                               self.frequency.read, self.brightness.read,
+                               p0=[self.flopping_freq])
+        perr = np.sqrt(np.diag(pcov))
+        if perr < 0.1:
+            self.flopping_freq = float(popt)
diff --git a/examples/pdb.pyon b/examples/pdb.pyon
index 0967ef424..4ca76a468 100644
--- a/examples/pdb.pyon
+++ b/examples/pdb.pyon
@@ -1 +1 @@
-{}
+{"flopping_freq": 1500.0294421161527}