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 = '../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 = 80 # from t
  24 tto = 100 # to t
  25 acangle = 58 # angle of ac drive
  26 dcangle = 0 # angle of dc bias
  27 sper = 4 # spacial periodicity
  28 pad = 0.01 # small padding so the ticks get drawn correctly...
  29 msize = 20
  30 res = (1440/80,900/80) # default dpi is 80
  31 cmap = cm.gray_r # colormap
  32 shader = False # applies ultra cool shadows to the map!
  33 show_plot = False
  34 xlabel = r'$x$'
  35 ylabel = r'$y$'
  36 cblabel = r'$\cos(x)\cos(y)+\cos(x)+\cos(y)$'
  37 fs_ticks = 16
  38 fs_labels = 24
  39 cof_only = False
  40 has_dc = False
  41 test = False
  42 title = 'F = 0.1, alpha = 83, l = 3.0, a = 1.5, eta1 = 1.0, eta2 = 1.5, theta = 0'
  43
  44 ################################################################################
  45
  46 print "reading data..."
  47
  48 data_cof= mlab.csv2rec(FCOF, delimiter='\t', comments='#')
  49 data_p1 = mlab.csv2rec(FP1, delimiter='\t', comments='#')
  50 data_p2 = mlab.csv2rec(FP2, delimiter='\t', comments='#')
  51
  52 t = data_cof.r
  53 cofx = data_cof.x
  54 cofy = data_cof.y
  55 p1x = data_p1.x
  56 p1y = data_p1.y
  57 p2x = data_p2.x
  58 p2y = data_p2.y
  59 pxdiff = data_p2.x - data_p1.x
  60 pydiff = data_p2.y - data_p1.y
  61
  62 print "done."
  63
  64 tstep = t[10] - t[9]
  65
  66 fromt = tfrom/tstep
  67 tot = tto/tstep
  68
  69 every = evry*tstep
  70 trunc_cofx = cofx[fromt*tau:tot*tau:evry].tolist()
  71 trunc_cofy = cofy[fromt*tau:tot*tau:evry].tolist()
  72 trunc_p1x = p1x[fromt*tau:tot*tau:evry].tolist()
  73 trunc_p2x = p2x[fromt*tau:tot*tau:evry].tolist()
  74 trunc_p2y = p2y[fromt*tau:tot*tau:evry].tolist()
  75 trunc_p1y = p1y[fromt*tau:tot*tau:evry].tolist()
  76 trunc_pxdiff = pxdiff[fromt*tau:tot*tau:evry].tolist()
  77 trunc_pydiff = pydiff[fromt*tau:tot*tau:evry].tolist()
  78 trunc_t = t[fromt*tau:tot*tau:evry].tolist()
  79
  80 trunc_p1x_eff = []
  81 trunc_p2x_eff = []
  82 trunc_p1y_eff = []
  83 trunc_p2y_eff = []
  84 trunc_pxdiff_eff = []
  85 trunc_pydiff_eff = []
  86
  87 for i in range(len(trunc_cofx)):
  88         assert len(trunc_cofx) == len(trunc_p1x)
  89         print i, len(trunc_p1x_eff), len(trunc_cofx)
  90         diff = abs(trunc_cofx[i]-trunc_cofx[i-1])
  91         if (diff < 1) and (i != 0) and i:
  92                 trunc_p1x[i] = trunc_p1x[i-1]
  93                 trunc_p2x[i] = trunc_p2x[i-1]
  94                 trunc_p1y[i] = trunc_p1y[i-1]
  95                 trunc_p2y[i] = trunc_p2y[i-1]
  96                 trunc_pxdiff[i] = trunc_pxdiff[i-1]
  97                 trunc_pydiff[i] = trunc_pydiff[i-1]
  98         else:
  99                         print diff
 100
 101 #print trunc_p1x_eff
 102
 103 xmax = np.ceil(max(trunc_cofx)/L)
 104 xmin = np.ceil(min(trunc_cofx)/L)
 105 ymax = np.ceil(max(trunc_cofy)/L)
 106 ymin = np.ceil(min(trunc_cofy)/L)
 107
 108 xdiff = np.ceil(abs(max(trunc_cofx)-min(trunc_cofx))/L)
 109 ydiff = np.ceil(abs(max(trunc_cofy)-min(trunc_cofy))/L)
 110
 111 xcorr = xdiff-sper
 112 ycorr = ydiff-sper
 113
 114 print xdiff, ydiff
 115 print xdiff-(xdiff-sper)
 116
 117 if xmax > 0:
 118         limx_min = (xmin+2)*L
 119         limx_max = (xmin+2+sper)*L
 120 if xmax < 0:
 121         limx_min = (xmax-1-sper)*L
 122         limx_max = (xmax-1)*L
 123 if ymax > 0:
 124         limy_min = (ymin-2)*L
 125         limy_max = (ymin-2+sper)*L
 126 if ymax < 0:
 127         limy_min = (ymax-sper)*L
 128         limy_max = (ymax)*L
 129
 130 print (limx_max-limx_min)/L
 131
 132 print "x:","(",limx_min,",",limx_max,")"
 133 print "y:","(",limy_min,",",limy_max,")"
 134
 135 xticks = np.arange(limx_min,limx_max+pad,L)
 136 yticks = np.arange(limy_min,limy_max+pad,L)
 137 xlim = (limx_min-pad,limx_max+pad)
 138 ylim = (limy_min-pad,limy_max+pad)
 139 tlabel = ("-2L","-L","0","L","2L")
 140
 141 print "generating basemap..."
 142 delta = 0.01
 143 x = np.arange(limx_min,limx_max, delta)
 144 y = np.arange(limy_min,limy_max, delta)
 145 X,Y = np.meshgrid(x, y)
 146 Z = np.cos(X)*np.cos(Y)+np.cos(X)+np.cos(Y)
 147 print "done."
 148
 149 if test is True:
 150         raise
 151
 152 fig = plt.figure(figsize=res)
 153
 154 ax = fig.add_subplot(111)
 155
 156 if (shader is True):
 157     print "applying lightsource..."
 158     ls = col.LightSource(azdeg=0,altdeg=90)
 159     rgb = ls.shade(Z,cmap)
 160     print "done."
 161     print "plotting basemap..."
 162     pot = plt.imshow(rgb, extent=[limx_min,limx_max,limy_min,limy_max])
 163     print "done."
 164 else:
 165     print "plotting basemap..."
 166     pot = plt.pcolormesh(X,Y,Z, cmap=cmap)
 167     print "done."
 168
 169 cb = plt.colorbar(pot)
 170 cb.set_label(cblabel)
 171
 172 print "plotting dimer onto the basemap..."
 173 if (cof_only is False):
 174 #       cof = ax.scatter(trunc_cofx, trunc_cofy, c='yellow', marker='o', lw=0, s=msize/4, label='cof')
 175         x1 = ax.scatter(trunc_p1x, trunc_p1y, c='black', marker='o', lw=0, s=msize, label='x1')
 176         x2 = ax.scatter(trunc_p2x, trunc_p2y, c='black', marker='o', lw=0, s=msize, label='x2')
 177
 178         k=0
 179         for i in range(len(trunc_p1x)):
 180                 if trunc_t[i] >= k*every:
 181                         #print "plotting because",t[i],"=",k,"*",every,"(",k*every,")"
 182                         bond = ax.arrow(trunc_p1x[i], trunc_p1y[i],
 183                                         trunc_pxdiff[i],
 184                                         trunc_pydiff[i], color='black',
 185                                         lw=0.1, label='bond',
 186                                         alpha=0.5)
 187                         k = k+1
 188                 elif trunc_t[i] < k*every:
 189                         #print "not plotting because",t[i],"!=",k,"*",every,"(",k*every,")"
 190                         continue
 191 elif (cof_only is True):
 192         cof = ax.scatter(trunc_cofx, trunc_cofy, c='red', marker='o',
 193                         lw=0, s=msize/4, label='cof')
 194 else:
 195         print "fatal error!"
 196
 197
 198 print "done."
 199
 200 ax.arrow(limx_min,limy_min, 10*np.cos(acangle*np.pi/180),
 201                 10*np.sin(acangle*np.pi/180), lw=2.5, color='blue',
 202                 head_width=0.15, shape='full')
 203
 204 if (has_dc is True):
 205         ax.arrow(limx_min,limy_min, 10*np.cos(dcangle*np.pi/180),
 206                         10*np.sin(dcangle*np.pi/180), lw=2.5, color='red',
 207                         head_width=0.15, shape='full')
 208
 209 ax.set_xlabel(xlabel, size=fs_labels)
 210 ax.set_ylabel(ylabel, size=fs_labels)
 211 ax.set_xticks(xticks)
 212 ax.set_yticks(yticks)
 213 ax.set_xlim(limx_min,limx_max)
 214 ax.set_ylim(limy_min,limy_max)
 215 ax.set_xticklabels(tlabel, size=fs_ticks)
 216 ax.set_yticklabels(tlabel, size=fs_ticks)
 217
 218 if (cof_only is False):
 219         legend = fig.legend((x1, x2, bond), (r'$\vec{r}_1$',
 220                 r'$\vec{r}_2$', r'$|\vec{r}_1-\vec{r}_2|$'), shadow=True,
 221                 fancybox=True)
 222 """
 223 elif (cof_only is True):
 224         legend = fig.legend((cof),(r'$\vec{r}_s$'), shadow=True, fancybox=True)
 225 else:
 226         print "fatal error!"
 227
 228 #legend.get_frame().set_alpha(0.75)
 229 """
 230
 231 #plt.title(title)
 232
 233 print "saving to file.."
 234 plt.savefig(FOUT)
 235 print "done."
 236
 237 if (show_plot is True):
 238         print "showing plot..."
 239         plt.show()
 240         print "done."
 241