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'$t$'
  19 ylabel1 = r'$<x^2>$'
  20 xlabel2 = r'$r$'
  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 len(tra_x[0])
  57 print len(x)
  58 print t
  59
  60 fig = plt.figure(figsize=res)
  61 ax1 = fig.add_subplot(121)
  62 ax2 = fig.add_subplot(122)
  63
  64
  65 D_x =[ [] for DUMMYVAR in range(len(unique(r))) ]
  66 D_y =[ [] for DUMMYVAR in range(len(unique(r))) ]
  67
  68 #D_x = []
  69 #D_y = []
  70
  71 #dmin = min(unique(D))
  72 #dmax = max(unique(D))
  73
  74
  75
  76 for i in range(len(unique(r))):
  77         for j in range(len(tra_x[i])):
  78         #       if (j > 0):
  79                 tmp = tra_x[i][j]/(2*j*tau)
  80                 D_x[i].append(tmp)
  81
  82 for i in range(len(unique(r))):
  83         for j in range(len(tra_y[i])):
  84         #       if (j > 0):
  85                 tmp = tra_y[i][j]/(2*j*tau)
  86                 D_y[i].append(tmp)
  87
  88 print len(D_x)
  89 print len(D_y)
  90
  91 print t
  92 print D_x[i]
  93
  94 for i in range(len(unique(r))):
  95         ax1.plot(t,tra_x[i])
  96
  97 for i in range(len(unique(r))):
  98         ax2.plot(t,D_x[i])
  99 #       ax2.plot(t,D_y[i])
 100
 101 ax1.set_xlim(0,len(tra_x[0]))
 102 #ax1.set_ylim(x.min(),x.max())
 103 ax1.set_xlabel(xlabel1)
 104 ax1.set_ylabel(ylabel1)
 105 ax2.set_xlim(0,len(tra_x[0]))
 106 #ax2.set_ylim(D.min(),D.max())
 107 ax2.set_xlabel(xlabel2)
 108 ax2.set_ylabel(ylabel2)
 109
 110 plt.savefig(FOUT)
 111 plt.show()