pinv2.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 = '../latex2/images/dimer_sympartner.png'
  13 FCOF = '../rk_out.dat'
  14 FP1 = '../rk_point1.dat'
  15 FP2 = '../rk_point2.dat'
  16 FVEL = "../rk_velo.dat"
  17
  18 L = 2*np.pi
  19 w = 0.1
  20 tau = L/w
  21 evry = 1 # plot periodicity, don't set < 10 (bug)
  22 offx = -98 # offset of (x,y) in L
  23 offy = 94
  24 tfrom = 80 # from t
  25 tto = 100 # to t
  26 acangle = 58 # angle of ac drive
  27 dcangle = 0 # angle of dc bias
  28 sperx = 4 # spacial periodicity
  29 spery = 3
  30 pad = 0.01 # small padding so the ticks get drawn correctly...
  31 msize = 40
  32 res = (1440/80,900/80) # default dpi is 80
  33 cmap = cm.gray_r # colormap
  34 shader = False # applies ultra cool shadows to the map!
  35 show_plot = False
  36 xlabel = r'$x$'
  37 ylabel = r'$y$'
  38 cblabel = r'$\cos(x)\cos(y)+\cos(x)+\cos(y)$'
  39 fs_ticks = 16
  40 fs_labels = 24
  41 cof_only = False
  42 has_dc = False
  43 test = False
  44 show_cb = False
  45 show_legend = False
  46 show_annotations = False
  47 thresh = 1
  48 title = 'F = 0.1, alpha = 83, l = 3.0, a = 1.5, eta1 = 1.0, eta2 = 1.5, theta = 0'
  49
  50 ################################################################################
  51
  52 print "reading data..."
  53
  54 data_cof= mlab.csv2rec(FCOF, delimiter='\t', comments='#')
  55 data_p1 = mlab.csv2rec(FP1, delimiter='\t', comments='#')
  56 data_p2 = mlab.csv2rec(FP2, delimiter='\t', comments='#')
  57 data_v = mlab.csv2rec(FVEL, delimiter='\t', comments='#')
  58
  59 t = data_cof.r
  60 cofx = data_cof.x
  61 cofy = data_cof.y
  62 p1x = data_p1.x
  63 p1y = data_p1.y
  64 p2x = data_p2.x
  65 p2y = data_p2.y
  66 pxdiff = data_p2.x - data_p1.x
  67 pydiff = data_p2.y - data_p1.y
  68 vx = data_v.x
  69 vy = data_v.y
  70
  71 print "done."
  72
  73 tstep = t[10] - t[9]
  74
  75 fromt = tfrom/tstep
  76 tot = tto/tstep
  77
  78 tcx = [0.59,0.65,0.65,0.57,0.52,0.39]
  79 tcy = [0.4,0.6,0.83,0.6,0.51,0.4]
  80 tl = ["1","2","3","4","5","6"]
  81
  82 every = evry*tstep
  83 trunc_cofx = cofx[fromt*tau:tot*tau:evry].tolist()
  84 trunc_cofy = cofy[fromt*tau:tot*tau:evry].tolist()
  85 trunc_p1x = p1x[fromt*tau:tot*tau:evry].tolist()
  86 trunc_p2x = p2x[fromt*tau:tot*tau:evry].tolist()
  87 trunc_p2y = p2y[fromt*tau:tot*tau:evry].tolist()
  88 trunc_p1y = p1y[fromt*tau:tot*tau:evry].tolist()
  89 trunc_pxdiff = pxdiff[fromt*tau:tot*tau:evry].tolist()
  90 trunc_pydiff = pydiff[fromt*tau:tot*tau:evry].tolist()
  91 trunc_t = t[fromt*tau:tot*tau:evry].tolist()
  92 trunc_vx = vx[fromt*tau:tot*tau:evry].tolist()
  93 trunc_vy = vy[fromt*tau:tot*tau:evry].tolist()
  94
  95 def uniquify(seq, idfun=None):
  96         # order preserving
  97         if idfun is None:
  98                 def idfun(x): return x
  99         seen = {}
 100         result = []
 101         for item in seq:
 102                 marker = idfun(item)
 103                 # in old Python versions:
 104                 # # if seen.has_key(marker)
 105                 # # but in new ones:
 106                 if marker in seen: continue
 107                 seen[marker] = 1
 108                 result.append(item)
 109         return result
 110
 111 trunc_p1x_eff = []
 112 trunc_p2x_eff = []
 113 trunc_p1y_eff = []
 114 trunc_p2y_eff = []
 115 trunc_pxdiff_eff = []
 116 trunc_pydiff_eff = []
 117 trunc_vx_eff = []
 118 trunc_vy_eff = []
 119
 120 for i in range(len(trunc_cofx)):
 121         assert len(trunc_cofx) == len(trunc_p1x)
 122         #print i, len(trunc_p1x_eff), len(trunc_cofx)
 123         vel = np.sqrt(pow(trunc_vx[i],2)+pow(trunc_vy[i],2))
 124         p1 = abs(np.sqrt(pow(trunc_p1x[i],2)+pow(trunc_p1y[i],2))-np.sqrt(pow(trunc_p1x[i-1],2)+pow(trunc_p1y[i-1],2)))
 125         p2 = abs(np.sqrt(pow(trunc_p2x[i],2)+pow(trunc_p2y[i],2))-np.sqrt(pow(trunc_p2x[i-1],2)+pow(trunc_p2y[i-1],2)))
 126         if (i != 0):
 127                 if (p1 < thresh) and (p2 < thresh):
 128                         trunc_cofx[i] = trunc_cofx[i-1]
 129                         trunc_cofy[i] = trunc_cofy[i-1]
 130                         trunc_p1x[i] = trunc_p1x[i-1]
 131                         trunc_p2x[i] = trunc_p2x[i-1]
 132                         trunc_p1y[i] = trunc_p1y[i-1]
 133                         trunc_p2y[i] = trunc_p2y[i-1]
 134                         trunc_pxdiff[i] = trunc_pxdiff[i-1]
 135                         trunc_pydiff[i] = trunc_pydiff[i-1]
 136                         trunc_vx[i] = trunc_vx[i-1]
 137                         trunc_vy[i] = trunc_vy[i-1]
 138         else:
 139                 print "velocity:",vel, "dp1:", p1, "dp2:", p2
 140
 141 set_p1x = trunc_p1x#uniquify(trunc_p1x)
 142 set_p1y = trunc_p1y#uniquify(trunc_p1y)
 143 set_p2x = trunc_p2x#uniquify(trunc_p2x)
 144 set_p2y = trunc_p2y#uniquify(trunc_p2y)
 145 set_pxdiff = trunc_pxdiff#uniquify(trunc_pxdiff)
 146 set_pydiff = trunc_pydiff#uniquify(trunc_pydiff)
 147 set_vx = trunc_vx#uniquify(trunc_vx)
 148 set_vy = trunc_vy#uniquify(trunc_vy)
 149 set_cofx = trunc_cofx
 150 set_cofy = trunc_cofy
 151
 152
 153 print len(set_cofx), len(set_cofy)
 154
 155 xmax = np.ceil(max(trunc_cofx)/L)
 156 xmin = np.ceil(min(trunc_cofx)/L)
 157 ymax = np.ceil(max(trunc_cofy)/L)
 158 ymin = np.ceil(min(trunc_cofy)/L)
 159
 160 xdiff = np.ceil(abs(max(trunc_cofx)-min(trunc_cofx))/L)
 161 ydiff = np.ceil(abs(max(trunc_cofy)-min(trunc_cofy))/L)
 162
 163 xcorr = xdiff-sperx
 164 ycorr = ydiff-spery
 165
 166 #print xdiff, ydiff
 167 #print xdiff-(xdiff-sper)
 168
 169 if xmax > 0:
 170         limx_min = (xmin+2)*L
 171         limx_max = (xmin+2+sperx)*L
 172 if xmax < 0:
 173         limx_min = (xmax-2-sperx)*L
 174         limx_max = (xmax-2)*L
 175 if ymax > 0:
 176         limy_min = (ymin-2)*L
 177         limy_max = (ymin-2+spery)*L
 178 if ymax < 0:
 179         limy_min = (ymax-spery)*L
 180         limy_max = (ymax)*L
 181
 182 #print (limx_max-limx_min)/L
 183
 184 print "x:","(",limx_min,",",limx_max,")"
 185 print "y:","(",limy_min,",",limy_max,")"
 186
 187 xticks = np.arange(limx_min,limx_max+pad,L)
 188 yticks = np.arange(limy_min,limy_max+pad,L)
 189 xlim = (limx_min-pad,limx_max+pad)
 190 ylim = (limy_min-pad,limy_max+pad)
 191 tlabel = ("-2L","-L","0","L","2L")
 192
 193 print "generating basemap..."
 194 delta = 0.01
 195 x = np.arange(limx_min,limx_max, delta)
 196 y = np.arange(limy_min,limy_max, delta)
 197 X,Y = np.meshgrid(x, y)
 198 Z = np.cos(X)*np.cos(Y)+np.cos(X)+np.cos(Y)
 199 print "done."
 200
 201 if test is True:
 202         raise
 203
 204 fig = plt.figure() #(figsize=res)
 205
 206 ax = fig.add_subplot(111)
 207
 208 if (shader is True):
 209     print "applying lightsource..."
 210     ls = col.LightSource(azdeg=0,altdeg=90)
 211     rgb = ls.shade(Z,cmap)
 212     print "done."
 213     print "plotting basemap..."
 214     pot = plt.imshow(rgb, extent=[limx_min,limx_max,limy_min,limy_max])
 215     print "done."
 216 else:
 217     print "plotting basemap..."
 218     pot = plt.pcolormesh(X,Y,Z, cmap=cmap)
 219     print "done."
 220
 221 if (show_cb is True):
 222         cb = plt.colorbar(pot)
 223         cb.set_label(cblabel)
 224
 225 print "plotting dimer onto the basemap..."
 226 if (cof_only is False):
 227         cof = ax.plot(trunc_cofx, trunc_cofy, c='green', marker='.',
 228                         lw=1, label='cof', alpha=0.5)
 229         x1 = ax.scatter(set_p1x, set_p1y, c='black', marker='o', lw=0, s=msize, label='x1')
 230         x2 = ax.scatter(set_p2x, set_p2y, c='grey', marker='o', lw=0, s=msize, label='x2')
 231
 232         k=0
 233         for i in range(len(set_p1x)):
 234                 if trunc_t[i] >= k*every:
 235                         #print "plotting because",t[i],"=",k,"*",every,"(",k*every,")"
 236                         bond = ax.arrow(set_p1x[i], set_p1y[i],
 237                                         set_pxdiff[i],
 238                                         set_pydiff[i], color='black',
 239                                         lw=1, label='bond',
 240                                         alpha=0.5)
 241                         velo = ax.arrow(set_cofx[i], set_cofy[i],
 242                                         2*set_vx[i], 2*set_vy[i],
 243                                         color='red', lw=1,
 244                                         head_width=0.15, shape='full',
 245                                         label='velo')
 246                         k = k+1
 247                 elif trunc_t[i] < k*every:
 248                         #print "not plotting because",t[i],"!=",k,"*",every,"(",k*every,")"
 249                         continue
 250 elif (cof_only is True):
 251         cof = ax.scatter(trunc_cofx, trunc_cofy, c='red', marker='o',
 252                         lw=0, s=msize/4, label='cof')
 253 else:
 254         print "fatal error!"
 255
 256
 257 print "done."
 258
 259 ax.arrow(limx_min+L/2,limy_min+L/2, 2*np.cos(acangle*np.pi/180),
 260                 2*np.sin(acangle*np.pi/180), lw=2.5, color='blue',
 261                 head_width=0.15, shape='full')
 262 ax.arrow(limx_min+L/2,limy_min+L/2, -2*np.cos(acangle*np.pi/180),
 263                 -2*np.sin(acangle*np.pi/180), lw=2.5, color='blue',
 264                 head_width=0.15, shape='full')
 265 if (has_dc is True):
 266         ax.arrow(limx_min+L/2,limy_min+L/2, 2*np.cos(dcangle*np.pi/180),
 267                         2*np.sin(dcangle*np.pi/180), lw=2.5, color='red',
 268                         head_width=0.15, shape='full')
 269
 270 if (show_annotations is True):
 271         for i in range(len(tcx)):
 272                 fig.text(tcx[i],tcy[i],tl[i],bbox=dict(boxstyle="round",
 273                         fc=(1.0, 0.7, 0.7), ec="none"))
 274
 275
 276 ax.set_xlabel(xlabel, size=fs_labels)
 277 ax.set_ylabel(ylabel, size=fs_labels)
 278 ax.set_xticks(xticks)
 279 ax.set_yticks(yticks)
 280 ax.set_xlim(limx_min,limx_max)
 281 ax.set_ylim(limy_min,limy_max)
 282 ax.set_xticklabels(tlabel, size=fs_ticks)
 283 ax.set_yticklabels(tlabel, size=fs_ticks)
 284
 285 if (show_legend is True):
 286         legend = fig.legend((x1, x2, bond), (r'$\vec{r}_1$',
 287                 r'$\vec{r}_2$', r'$|\vec{r}_1-\vec{r}_2|$'), shadow=True,
 288                 fancybox=True)
 289 """
 290 elif (cof_only is True):
 291         legend = fig.legend((cof),(r'$\vec{r}_s$'), shadow=True, fancybox=True)
 292 else:
 293         print "fatal error!"
 294
 295 #legend.get_frame().set_alpha(0.75)
 296 """
 297
 298 #plt.title(title)
 299
 300 print "saving to file.."
 301 plt.savefig(FOUT)
 302 print "done."
 303
 304 if (show_plot is True):
 305         print "showing plot..."
 306         plt.show()
 307         print "done."
 308
 309