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