Blame view

ringdown.py 1.12 KB
84d931fb3   bmarechal   .
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
  import matplotlib.pyplot as plt
  import numpy as np
  import csv, glob, sys
  from scipy.optimize import curve_fit
  
  list_files = (glob.glob(sys.argv[1]))
  
  data = []
  
  for f in list_files:
      data_iter = csv.reader(open(f, 'r'), delimiter = ',', quotechar = '"')
      for i in range(2):
          data_iter.next()
      #temp_data = [value for value in data_iter]      #if ends with a number
      temp_data = [value[:-1] for value in data_iter]#if ends with a comma
      data.extend(temp_data)
  
  data = np.asarray(data, dtype = float)
  
  del(temp_data, list_files, value, f)
  
  plt.subplot(111)
  plt.clf()
  
  plt.plot(data[data[:,0]>=0,0], data[data[:,0]>=0,1], label ='mes')
  
  def func(t, tau, A, w , a, b):
      return A * np.exp(-t/tau) * np.sin((w+a*t)*t) + b
  
  popt, pcov = curve_fit(func, data[data[:,0]>=0,0], data[data[:,0]>=0,1], p0 = [1e-5, 1e-2, 1e5, 1e10, 1e-4], maxfev=10000)
  yfit = func(data[data[:,0]>=0,0], *popt)
  
  tau = float(popt[0])
  Q = np.pi*299792458*tau/(2.*140e-3)
  
  #print(np.sqrt(np.diag(pcov)))
  print(tau, Q)
  
  plt.plot(data[data[:,0]>=0,0], yfit, label ='fit')
  
  plt.xlabel('t')
  plt.ylabel('Intensity')
  plt.grid(which='both')
  plt.legend()
  plt.show()