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