detdiff.py

   1 import matplotlib as mpl
   2 import matplotlib.mlab as mlab
   3 import matplotlib.pyplot as plt
   4 import matplotlib.axes as axe
   5 import matplotlib.collections as collections
   6 import numpy as np
   7 import random
   8 from pylab import *
   9 from matplotlib.pyplot import *
  10
  11 FIN = '../rk_variance.dat'
  12 FOUT = '../diffusion.png'
  13
  14 L = 2*np.pi;
  15 tau = L/0.1;
  16 samples = 3
  17
  18 xlabel1 = r'$k\,[\tau]$'
  19 ylabel1 = r'$\langle {x(k\tau)}^2 - \langle x(k\tau) \rangle^2\rangle$'
  20 xlabel2 = r'$k\,[\tau]$'
  21 ylabel2 = r'$D$'
  22 res = (1440/80,900/80) # default dpi is 80
  23
  24 def get_trajectory(r,x):
  25         print "---------------------------------------------------------"
  26         uniquer = unique(r)
  27         ret_x = [ [] for DUMMYVAR in range(len(uniquer)) ]
  28         ret_y = [ [] for DUMMYVAR in range(len(uniquer)) ]
  29         for j in range(len(uniquer)):
  30                 k = 0
  31                 for i in range(len(x)):
  32                         if (r[i] == uniquer[j]):
  33                                 k = k+1
  34                                 ret_x[j].append(x[i])
  35                                 #ret_y[j].append(y[i])
  36                 print uniquer[j],":",len(ret_x[j]),k
  37         print "---------------------------------------------------------"
  38         #return (ret_x, ret_y)
  39         return (ret_x)
  40
  41 print "reading data..."
  42
  43 data  = mlab.csv2rec(FIN, comments='#', delimiter='\t')
  44 r = data.r
  45 x = data.x
  46 y = data.y
  47
  48 print "done."
  49
  50 print "extracting trajectories..."
  51 tra_x = get_trajectory(r,x)
  52 tra_y = get_trajectory(r,y)
  53 print "done."
  54
  55 t = range(len(tra_x[0]))
  56 print t
  57
  58 fig = plt.figure(figsize=res)
  59 ax1 = fig.add_subplot(121)
  60 ax2 = fig.add_subplot(122)
  61
  62
  63 D_x =[ [] for DUMMYVAR in range(len(unique(r))) ]
  64 D_y =[ [] for DUMMYVAR in range(len(unique(r))) ]
  65
  66 #D_x = []
  67 #D_y = []
  68
  69 #dmin = min(unique(D))
  70 #dmax = max(unique(D))
  71
  72
  73
  74 for i in range(len(unique(r))):
  75         for j in t:
  76                 tmp = tra_x[i][j]/(2*(j+1)*tau)
  77                 #print tmp,"=",tra_x[i][j],"/",2*(j+1)*tau
  78                 D_x[i].append(tmp)
  79
  80
  81 for i in range(len(unique(r))):
  82         for j in t:
  83                 tmp = tra_y[i][j]/(2*(j+1)*tau)
  84                 D_y[i].append(tmp)
  85
  86 #print tra_x
  87 #print D_x
  88
  89
  90 for i in range(len(unique(r))):
  91         ax1.plot(t,tra_x[i])
  92
  93 for i in range(len(unique(r))):
  94         ax2.plot(t,D_x[i])
  95 #       ax2.plot(t,D_y[i])
  96
  97 #ax1.set_xlim(0,len(tra_x[0]))
  98 ax1.set_ylim(0,14001)
  99 ax1.set_xlabel(xlabel1)
 100 ax1.set_ylabel(ylabel1)
 101 #ax2.set_xlim(0,len(tra_x[0]))
 102 ax2.set_ylim(0,2)
 103 ax2.set_xlabel(xlabel2)
 104 ax2.set_ylabel(ylabel2)
 105
 106 plt.savefig(FOUT)
 107 #plt.show()