pinv.py

   1 #!/usr/bin/python
   2 # -*- coding: iso-8859-1 -*-
   3
   4 import matplotlib as mpl
   5 import numpy as np
   6 import matplotlib.cm as cm
   7 import matplotlib.mlab as mlab
   8 import matplotlib.pyplot as plt
   9 import matplotlib.patches as patch
  10 import matplotlib.colors as col
  11
  12 FOUT = '../dimer_in_potential.png'
  13 FCOF = '../rk_out.dat'
  14 FP1 = '../rk_point1.dat'
  15 FP2 = '../rk_point2.dat'
  16
  17 L = 2*np.pi
  18 w = 0.1
  19 tau = L/w
  20 evry = 1 # plot periodicity, don't set < 10 (bug)
  21 offx = -98 # offset of (x,y) in L
  22 offy = 94
  23 tfrom = 90 # from t
  24 tto = 100 # to t
  25 driveangle = 41 # angle of ac drive
  26 sper = 4 # spacial periodicity
  27 pad = 0.01 # small padding so the ticks get drawn correctly...
  28 msize = 20
  29 res = (1440/80,900/80) # default dpi is 80
  30 cmap = cm.gray_r # colormap
  31 shader = False # applies ultra cool shadows to the map!
  32 show_plot = False
  33 xlabel = r'$x$'
  34 ylabel = r'$y$'
  35 cblabel = r'$\cos(x)\cos(y)+\cos(x)+\cos(y)$'
  36 fs_ticks = 14
  37 fs_labels = 24
  38 title = 'F = 0.1, alpha = 83, l = 3.0, a = 1.5, eta1 = 1.0, eta2 = 1.5, theta = 0'
  39
  40 ################################################################################
  41
  42 print "reading data..."
  43
  44 data_cof= mlab.csv2rec(FCOF, delimiter='\t', comments='#')
  45 data_p1 = mlab.csv2rec(FP1, delimiter='\t', comments='#')
  46 data_p2 = mlab.csv2rec(FP2, delimiter='\t', comments='#')
  47
  48 t = data_cof.r
  49 cofx = data_cof.x
  50 cofy = data_cof.y
  51 p1x = data_p1.x
  52 p1y = data_p1.y
  53 p2x = data_p2.x
  54 p2y = data_p2.y
  55 pxdiff = data_p2.x - data_p1.x
  56 pydiff = data_p2.y - data_p1.y
  57
  58 print "done."
  59
  60 tstep = t[10] - t[9]
  61
  62 fromt = tfrom/tstep
  63 tot = tto/tstep
  64
  65 every = evry*tstep
  66 trunc_cofx = cofx[fromt*tau:tot*tau:evry]
  67 trunc_cofy = cofy[fromt*tau:tot*tau:evry]
  68 trunc_p1x = p1x[fromt*tau:tot*tau:evry]
  69 trunc_p2x = p2x[fromt*tau:tot*tau:evry]
  70 trunc_p2y = p2y[fromt*tau:tot*tau:evry]
  71 trunc_p1y = p1y[fromt*tau:tot*tau:evry]
  72 trunc_pxdiff = pxdiff[fromt*tau:tot*tau:evry]
  73 trunc_pydiff = pydiff[fromt*tau:tot*tau:evry]
  74 trunc_t = t[fromt*tau:tot*tau:evry]
  75
  76 xmax = np.ceil(max(trunc_cofx)/L)
  77 xmin = np.ceil(min(trunc_cofx)/L)
  78 ymax = np.ceil(max(trunc_cofy)/L)
  79 ymin = np.ceil(min(trunc_cofy)/L)
  80
  81 print xmin*L, xmax*L
  82 xdiff = np.ceil(abs(max(trunc_cofx)-min(trunc_cofx))/L)
  83 ydiff = np.ceil(abs(max(trunc_cofy)-min(trunc_cofy))/L)
  84
  85 xcorr = xdiff-sper
  86 ycorr = ydiff-sper
  87
  88 print xdiff, ydiff
  89 print xdiff-(xdiff-sper)
  90
  91 if xmax > 0:
  92         limx_max = xmax*L
  93         limx_min = (xmax-sper)*L
  94 if xmax < 0:
  95         limx_min = xmin*L
  96         limx_max = (xmin+sper)*L
  97 if ymax > 0:
  98         limy_max = ymax*L
  99         limy_min = (ymax-sper)*L
 100 if ymax < 0:
 101         limy_min = ymin*L
 102         limy_max = (ymin+sper)*L
 103
 104 print (limx_max-limx_min)/L
 105
 106 print "x:","(",limx_min,",",limx_max,")"
 107 print "y:","(",limy_min,",",limy_max,")"
 108
 109 xticks = np.arange(limx_min,limx_max+pad,L)
 110 yticks = np.arange(limy_min,limy_max+pad,L)
 111 xlim = (limx_min-pad,limx_max+pad)
 112 ylim = (limy_min-pad,limy_max+pad)
 113 tlabel = ("-2L","-L","0","L","2L")
 114
 115 print "generating basemap..."
 116 delta = 0.01
 117 x = np.arange(limx_min,limx_max, delta)
 118 y = np.arange(limy_min,limy_max, delta)
 119 X,Y = np.meshgrid(x, y)
 120 Z = np.cos(X)*np.cos(Y)+np.cos(X)+np.cos(Y)
 121 print "done."
 122
 123 fig = plt.figure(figsize=res)
 124
 125 ax = fig.add_subplot(111)
 126
 127 if (shader is True):
 128     print "applying lightsource..."
 129     ls = col.LightSource(azdeg=0,altdeg=90)
 130     rgb = ls.shade(Z,cmap)
 131     print "done."
 132     print "plotting basemap..."
 133     pot = plt.imshow(rgb, extent=[limx_min,limx_max,limy_min,limy_max])
 134     print "done."
 135 else:
 136     print "plotting basemap..."
 137     pot = plt.pcolormesh(X,Y,Z, cmap=cmap)
 138     print "done."
 139
 140 cb = plt.colorbar(pot)
 141 cb.set_label(cblabel)
 142
 143 print "plotting dimer onto the basemap..."
 144 cof = ax.scatter(trunc_cofx, trunc_cofy, c='yellow', marker='o', lw=0, s=msize/4, label='cof')
 145 x1 = ax.scatter(trunc_p1x, trunc_p1y, c='blue', marker='o', lw=0, s=msize, label='x1')
 146 x2 = ax.scatter(trunc_p2x, trunc_p2y, c='purple', marker='o', lw=0, s=msize, label='x2')
 147 #bond = ax.arrow(trunc_p1x, trunc_p1y, trunc_pxdiff, trunc_pydiff, color='r', lw=0.1, label='bond')
 148
 149
 150 k=0
 151 for i in range(len(trunc_t)):
 152         #bond = ax.arrow(p1x[i], p1y[i], pxdiff[i], pydiff[i], color='r', lw=0.1, label='bond')
 153         if trunc_t[i] >= k*every:
 154                 #print "plotting because",t[i],"=",k,"*",every,"(",k*every,")"
 155                 bond = ax.arrow(trunc_p1x[i], trunc_p1y[i],
 156                                 trunc_pxdiff[i], trunc_pydiff[i],
 157                                 color='r', lw=0.1, label='bond',
 158                                 alpha=0.5)
 159                 k = k+1
 160         elif trunc_t[i] < k*every:
 161                 #print "not plotting because",t[i],"!=",k,"*",every,"(",k*every,")"
 162                 continue
 163
 164 print "done."
 165
 166 ax.arrow(limx_min,limy_min, 10*np.cos(driveangle*np.pi/180),
 167                 10*np.sin(driveangle*np.pi/180), lw=5.0, color='g')
 168
 169 """
 170 n = 0
 171 for i in range(len(x)):
 172         if x[i] >= n*L:
 173                 s1 = ax.axhline(y=n/2)
 174                 n = n+1
 175         elif x[i] < n*L:
 176                 continue
 177 """
 178
 179
 180 ax.set_xlabel(xlabel, size=fs_labels)
 181 ax.set_ylabel(ylabel, size=fs_labels)
 182 ax.set_xticks(xticks)
 183 ax.set_yticks(yticks)
 184 ax.set_xlim(limx_min,limx_max)
 185 ax.set_ylim(limy_min,limy_max)
 186 ax.set_xticklabels(tlabel, size=fs_ticks)
 187 ax.set_yticklabels(tlabel, size=fs_ticks)
 188
 189 legend = fig.legend((cof, x1, x2, bond), (r'$\vec{r}_s$', r'$\vec{r}_1$',
 190         r'$\vec{r}_2$', r'$|\vec{r}_1-\vec{r}_2|$'), shadow=True,
 191         fancybox=True)
 192 legend.get_frame().set_alpha(0.75)
 193
 194 #plt.title(title)
 195
 196 print "saving to file.."
 197 plt.savefig(FOUT)
 198 print "done."
 199
 200 if (show_plot is True):
 201         print "showing plot..."
 202         plt.show()
 203         print "done."
 204