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 FVAR = '../data/rk_variance.dat'
  12 #FMEAN = '../rk_mean.dat'
  13 FOUT = '../latex/images/diffusion.png'
  14
  15 L = 2*np.pi;
  16 tau = L/0.1;
  17 samples = 3
  18
  19 ylabel1 = r'$k\,[\tau]$'
  20 #ylabel1 = r'$\langle {x(k\tau)}^2 - \langle x(k\tau) \rangle^2\rangle$'
  21 xlabel1 = r'$\langle x \rangle$'
  22 xlabel2 = r'$k\,[\tau]$'
  23 ylabel2 = r'$D$'
  24 res = (1024/80,768/80) # default dpi is 80
  25
  26 def get_trajectory(r,x):
  27         print "---------------------------------------------------------"
  28         uniquer = unique(r)
  29         ret_x = [ [] for DUMMYVAR in range(len(uniquer)) ]
  30         ret_y = [ [] for DUMMYVAR in range(len(uniquer)) ]
  31         for j in range(len(uniquer)):
  32                 k = 0
  33                 for i in range(len(x)):
  34                         if (r[i] == uniquer[j]):
  35                                 k = k+1
  36                                 ret_x[j].append(x[i])
  37                                 #ret_y[j].append(y[i])
  38                 print uniquer[j],":",len(ret_x[j]),k
  39         print "---------------------------------------------------------"
  40         #return (ret_x, ret_y)
  41         return (ret_x)
  42
  43 print "reading data..."
  44
  45 var  = mlab.csv2rec(FVAR, comments='#', delimiter='\t')
  46 #mean = mlab.csv2rec(FMEAN, comments='#', delimiter='\t')
  47
  48 r = var.r
  49 var_x = var.x
  50 var_y = var.y
  51 #mean_x = mean.x
  52 #mean_y = mean.y
  53
  54 print "done."
  55
  56 print "extracting trajectories..."
  57 vartra_x = get_trajectory(r,var_x)
  58 vartra_y = get_trajectory(r,var_y)
  59 #meantra_x = get_trajectory(r,mean_x)
  60 #meantra_y = get_trajectory(r,mean_y)
  61 print "done."
  62
  63 t = range(len(vartra_x[0]))
  64 print t
  65 unique_r = unique(r)
  66
  67 fig = plt.figure(figsize=res)
  68 #ax1 = fig.add_subplot(121)
  69 ax2 = fig.add_subplot(111)
  70
  71
  72 D_x =[ [] for DUMMYVAR in range(len(unique(r))) ]
  73 D_y =[ [] for DUMMYVAR in range(len(unique(r))) ]
  74
  75 #D_x = []
  76 #D_y = []
  77
  78 #dmin = min(unique(D))
  79 #dmax = max(unique(D))
  80
  81
  82
  83 for i in range(len(unique(r))):
  84         for j in t:
  85                 tmp = vartra_x[i][j]/(2*(j+1)*tau)
  86                 #print tmp,"=",tra_x[i][j],"/",2*(j+1)*tau
  87                 D_x[i].append(tmp)
  88
  89
  90 for i in range(len(unique(r))):
  91         for j in t:
  92                 tmp = vartra_y[i][j]/(2*(j+1)*tau)
  93                 D_y[i].append(tmp)
  94
  95 #print tra_x
  96 #print D_x
  97
  98
  99 #for i in range(len(unique(r))):
 100 #       ax1.plot(meantra_x[i],t)
 101
 102 for i in range(len(unique(r))):
 103         ax2.plot(t,D_x[i])
 104 #       ax2.plot(t,D_y[i])
 105
 106 for i in range(len(unique(r))):
 107                 #tmp_tra = meantra_x[i]
 108                 tmp_D = D_x[i]
 109                 #ax1.text(tmp_tra[max(t)],max(t)+3,unique_r[i],)
 110                 ax2.text(max(t)+3,tmp_D[max(t)],unique_r[i],)
 111
 112 #ax1.set_xlim(0,len(tra_x[0]))
 113 #ax1.set_ylim(0,14001)
 114 #ax1.set_xlabel(xlabel1)
 115 #ax1.set_ylabel(ylabel1)
 116 #ax2.set_xlim(0,len(tra_x[0]))
 117 ax2.set_ylim(0,2)
 118 ax2.set_xlabel(xlabel2)
 119 ax2.set_ylabel(ylabel2)
 120
 121 plt.savefig(FOUT)
 122 #plt.show()