From a17e005aa906d6c1961c8cc36acd21397d2bec5b Mon Sep 17 00:00:00 2001
From: Rainer Wittmaack <ningo@gmx.net>
Date: Fri, 2 Oct 2009 10:56:55 +0200
Subject: [PATCH] * dimer.py -- plots a default dimer * fixed calculation of
 the deflection angle * enhanced the deflection angle plot * minor tweaks to
 plot fonts

---
 bifandvel.py   |   6 +--
 calculator.cpp |  23 ++++++--
 defangle.py    |  22 +++++---
 dimer.py       | 164 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 pinv.py        |  14 +++--
 5 files changed, 209 insertions(+), 20 deletions(-)
 create mode 100644 dimer.py

diff --git a/bifandvel.py b/bifandvel.py
index 2dd1bac..ca5547b 100644
--- a/bifandvel.py
+++ b/bifandvel.py
@@ -7,8 +7,8 @@ import matplotlib.mlab as mlab
 import matplotlib.axes as axe
 import matplotlib.text as txt
 
-DATAX = '../rk_bifurcate.dat'
-DATAV = '../rk_bif_avg_velo.dat'
+DATAX = '/home/raw/rk_bifurcate.dat'
+DATAV = '/home/raw/rk_bif_avg_velo.dat'
 FOUT = '../test.png'
 
 # some macros
@@ -20,7 +20,7 @@ xpad = 0.002
 ypad = 0.1
 max = 0.875
 tsize = 24
-ticksize = 14
+ticksize = 16
 
 # default values
 xtks1 = []
diff --git a/calculator.cpp b/calculator.cpp
index 80dc412..d8a1e1e 100644
--- a/calculator.cpp
+++ b/calculator.cpp
@@ -523,8 +523,11 @@ void Calculator::bifparam(const int &v, const QString &param, const Doub &r,
             params[7] = params[3]*params[11]*params[5]; // l1;
             params[8] = params[3]*params[11]*params[6]; // l2;
         }
-        else if (param == "t")
+        else if (param == "t") {
             params[10] = fmod(r,360)*PI/180;
+            params[15] = cos(params[10])*params[2];
+            params[16] = sin(params[10])*params[2];
+        }
         else
             throw("Fatal error!");
         break;
@@ -587,7 +590,7 @@ void Calculator::calc2(const Int n)
     qDebug() << "System" << v;
 
     Doub x, hnew;
-    Doub trajectory_length, deflection;
+    Doub trajectory_length, deflection, velocity_angle, bias_angle;
     Int k;
 
     Doub xpi[samples][s];
@@ -815,13 +818,23 @@ void Calculator::calc2(const Int n)
                         case 4:
                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
                         // makes only sense with r = params[10]
-                        trajectory_length = sqrt(SQR(abs(yout[0]-x0))+SQR(abs(yout[1]-y0)));
-                        deflection = ABS(fmod(acos(abs(yout[0]-x0)/trajectory_length)-params[10],(2*PI)));
+                        trajectory_length = sqrt(SQR(yout[0]-x0)+SQR(yout[1])-y0);
+                        if (abs(yout[0]-x0) > trajectory_length) {
+                            velocity_angle = PI;
+                        } else {
+                            velocity_angle = (acos((yout[0]-x0)/trajectory_length));
+                        }
+                        bias_angle = params[10];
+                        deflection = (180/PI)*ABS(velocity_angle - bias_angle);
+                        qDebug() << velocity_angle << "("<<yout[0]-x0<<")"
+                                << "-" << bias_angle << "="
+                                << deflection
+                                << "[ length:"<<trajectory_length<<"]";
                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
                         rkAvgvelo_x.insert(r,avg[0]);
                         rkAvgvelo_y.insert(r,avg[1]);
                         rkAvgvelo_phi.insert(r,avg[2]);
-                        rkDefAngle.insert(r,(180/PI)*deflection);
+                        rkDefAngle.insert(r,deflection);
                         break;
                     } // end switch
                 } // samples loop
diff --git a/defangle.py b/defangle.py
index 04b7cd0..0ee3c4b 100644
--- a/defangle.py
+++ b/defangle.py
@@ -14,6 +14,9 @@ FOUT = "../defangle.png"
 
 xlabel = r'$h$'
 ylabel = r'$\delta y$'
+res = (1440/80,900/80)
+pad1 = 0.01
+pad2 = 0.01
 
 ######################################################################
 
@@ -34,26 +37,29 @@ def mean_angle(param,angle):
 		for j in range(len(angle)):
 			if (param[j] == unique_param[i]):
 				list_angles[i].append(angle[j])
-				print len(list_angles[i])
+				#print len(list_angles[i])
 		ret_avg.append(float(sum(list_angles[i])) / len(list_angles[i]))
 		#print len(ret_avg)
 	return ret_avg
 
 
-fig = plt.figure()
+fig = plt.figure(figsize=res)
 
 ax1 = fig.add_subplot(211)
 ax2 = fig.add_subplot(212)
 
 avg_angles = mean_angle(r,psi)
+avg_velo = mean_angle(r,vx)
 unique_r = unique(r)
 
-ax1.plot(unique_r, avg_angles, marker=',', ms=0.1, ls='')
-ax2.plot(r, vx, marker=',', ms=0.1, ls='')
+ax1.plot(unique_r, avg_angles, marker='.', color='black')
+ax2.plot(unique_r, avg_velo, marker='.', color='black')
 
-ax1.set_xlim(r.min(),r.max())
-ax1.set_ylim(psi.min(),psi.max())
+ax1.set_xlim(r.min()-pad1,r.max()+pad1)
+ax1.set_ylim(psi.min()-pad1,psi.max()+pad1)
+ax2.set_xlim(r.min()-pad2,r.max()+pad2)
+ax2.set_ylim(vx.min()-pad2,vx.max()+pad2)
 
-#plt.savefig(FOUT)
+plt.savefig(FOUT)
 
-plt.show()
+#plt.show()
diff --git a/dimer.py b/dimer.py
new file mode 100644
index 0000000..3143a50
--- /dev/null
+++ b/dimer.py
@@ -0,0 +1,164 @@
+#!/usr/bin/python
+# -*- coding: iso-8859-1 -*-
+
+import matplotlib as mpl
+mpl.rc('text', usetex = True)
+import pylab
+import numpy as np
+import matplotlib.cm as cm
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+import matplotlib.patches as patch
+import matplotlib.colors as col
+
+FOUT = "../dimer_intro.png"
+
+L = 2*np.pi
+w = 0.1
+tau = L/w
+angle1 = 30
+angle2 = 60
+transport_angle = 100
+phi = 0
+alpha = 135
+dimer_length = 0.4
+cofx = cofy = L/2
+pad = 0.01 # small padding so the ticks get drawn correctly...
+msize = 20
+hw = 0.15 # arrow head width
+res = (1440/80,900/80) # default dpi is 80
+cmap = cm.gray_r # colormap
+shader = False # applies ultra cool shadows to the map!
+show_plot = False
+xticks = [-2*L,-L,0,L,2*L]
+yticks = [-2*L,-L,0,L,2*L]
+tlabel = ["-2L","-L","0","L","2L"]
+xlabel = r'$x$'
+ylabel = r'$y$'
+cblabel = r'$\cos(x)\cos(y)+\cos(x)+\cos(y)$'
+fs_ticks = 14
+fs_labels = 24
+
+################################################################################
+
+PI = np.pi
+rad1 = angle1*PI/180
+rad2 = angle2*PI/180
+rad_transport = transport_angle*PI/180
+rad_phi = phi*PI/180
+rad_alpha = alpha*PI/180
+
+limx_min = -L/2
+limx_max = 1.5*L
+# quadratic
+limy_min = limx_min
+limy_max = limx_max
+
+print "generating basemap..."
+delta = 0.01
+x = np.arange(limx_min,limx_max, delta)
+y = np.arange(limy_min,limy_max, delta)
+X,Y = np.meshgrid(x, y)
+Z = np.cos(X)*np.cos(Y)+np.cos(X)+np.cos(Y)
+print "done."
+
+fig = plt.figure(figsize=res)
+
+ax = fig.add_subplot(111)
+
+if (shader is True):
+    print "applying lightsource..."
+    ls = col.LightSource(azdeg=0,altdeg=90)
+    rgb = ls.shade(Z,cmap)
+    print "done."
+    print "plotting basemap..."
+    pot = plt.imshow(rgb, extent=[limx_min,limx_max,limy_min,limy_max])
+    print "done."
+else:
+    print "plotting basemap..."
+    pot = plt.pcolormesh(X,Y,Z, cmap=cmap)
+    print "done."
+    
+cb = plt.colorbar(pot)
+cb.set_label(cblabel)
+
+print "plotting dimer onto the basemap..."
+
+def make_circle(start,stop,r):
+	t = np.arange(start,stop,0.01)
+	#t = t.reshape((len(t), 1))
+	ret_x = cofx + r * np.cos(t)
+	ret_y = cofy + r * np.sin(t)
+	return ret_x, ret_y
+
+
+varphi_plot_x, varphi_plot_y = make_circle(rad1, rad2, 2.0)
+psi_plot_x, psi_plot_y = make_circle(rad_phi, rad_transport, 0.8)
+
+ax.scatter(cofx, cofy, c='yellow', marker='o', lw=0, s=msize*4,
+label='cof')
+ax.arrow(cofx, cofy, np.cos(rad1)*dimer_length*L,
+		np.sin(rad1)*dimer_length*L, color='r', lw=0.1, alpha=0.5)
+ax.arrow(cofx, cofy, np.cos(rad2)*dimer_length*L,
+		np.sin(rad2)*dimer_length*L, color='r', lw=0.1, alpha=0.5)
+ax.arrow(cofx, cofy, -np.cos(rad1)*dimer_length*L,
+		-np.sin(rad1)*dimer_length*L, color='r', lw=0.1, alpha=0.5)
+ax.arrow(cofx, cofy, -np.cos(rad2)*dimer_length*L,
+		-np.sin(rad2)*dimer_length*L, color='r', lw=0.1, alpha=0.5)
+ax.scatter(cofx+np.cos(rad1)*dimer_length*L,
+		cofy+np.sin(rad1)*dimer_length*L, c='blue', marker='o',
+		lw=0, s=msize, label='x1')
+ax.scatter(cofx-np.cos(rad1)*dimer_length*L,
+		cofy-np.sin(rad1)*dimer_length*L, c='purple', marker='o',
+		lw=0, s=msize, label='x1')
+ax.scatter(cofx+np.cos(rad2)*dimer_length*L,
+		cofy+np.sin(rad2)*dimer_length*L, c='blue',
+		marker='o', lw=0, s=msize, label='x2')
+ax.scatter(cofx-np.cos(rad2)*dimer_length*L,
+		cofy-np.sin(rad2)*dimer_length*L, c='purple',
+		marker='o', lw=0, s=msize, label='x2')
+ax.plot(varphi_plot_x, varphi_plot_y, color='black')
+ax.arrow(cofx, cofy, np.cos(rad_phi)*0.4*L, np.sin(rad_phi)*0.4*L,
+		head_width=hw, shape='full', color='red')
+ax.arrow(cofx, cofy, np.cos(rad_transport)*0.4*L,
+		np.sin(rad_transport)*0.4*L, head_width=hw,
+		shape='full')
+ax.plot(psi_plot_x, psi_plot_y, color='black')
+ax.arrow(cofx, cofy, np.cos(rad_alpha)*0.4*L, np.sin(rad_alpha)*0.4*L,
+		head_width=hw, shape='full', color='blue')
+ax.arrow(cofx, cofy, -np.cos(rad_alpha)*0.4*L, -np.sin(rad_alpha)*0.4*L,
+		head_width=hw, shape='full', color='blue')
+
+plt.text(cofx+np.cos(rad_phi+0.025)*0.4*L,
+		cofy+np.sin(rad_phi+0.025)*0.4*L,
+		r'$F\vec{e}_{\phi}$',{'color' : 'k', 'fontsize' :
+			fs_labels})
+plt.text(cofx+np.cos(rad_transport-0.05)*0.4*L,
+		cofy+np.sin(rad_transport-0.05)*0.4*L, r'$\vec{v}$',
+		{'color': 'k', 'fontsize' : fs_labels})
+plt.text(0.5*L,0.65*L, r'$\psi$', {'color' : 'k', 'fontsize' : fs_labels})
+plt.text(0.65*L,0.65*L, r'$\varphi$', {'color' : 'k', 'fontsize' :
+	fs_labels})
+plt.text(cofx+np.cos(rad_alpha-0.025)*0.4*L, cofy+np.sin(rad_alpha-0.025)*0.4*L,
+		r'$A\vec{e}_{\alpha}$', {'color' : 'k', 'fontsize' :
+			fs_labels})
+
+ax.set_xlabel(xlabel, size=fs_labels)
+ax.set_ylabel(ylabel, size=fs_labels)
+ax.set_xticks(xticks)
+ax.set_yticks(yticks)
+ax.set_xlim(limx_min,limx_max)
+ax.set_ylim(limy_min,limy_max)
+ax.set_xticklabels(tlabel, size=fs_ticks)
+ax.set_yticklabels(tlabel, size=fs_ticks)
+
+
+print "saving to file.."
+plt.savefig(FOUT)
+print "done."
+
+if (show_plot is True):
+	print "showing plot..."
+	plt.show()
+	print "done."
+
diff --git a/pinv.py b/pinv.py
index ec419b8..a01d0f2 100644
--- a/pinv.py
+++ b/pinv.py
@@ -22,7 +22,8 @@ offx = -98 # offset of (x,y) in L
 offy = 94
 tfrom = 90 # from t
 tto = 100 # to t
-driveangle = 41 # angle of ac drive
+acangle = 58 # angle of ac drive
+dcangle = 0 # angle of dc bias
 sper = 4 # spacial periodicity
 pad = 0.01 # small padding so the ticks get drawn correctly...
 msize = 20
@@ -33,7 +34,7 @@ show_plot = False
 xlabel = r'$x$'
 ylabel = r'$y$'
 cblabel = r'$\cos(x)\cos(y)+\cos(x)+\cos(y)$'
-fs_ticks = 14
+fs_ticks = 16
 fs_labels = 24
 title = 'F = 0.1, alpha = 83, l = 3.0, a = 1.5, eta1 = 1.0, eta2 = 1.5, theta = 0'
 
@@ -163,8 +164,13 @@ for i in range(len(trunc_t)):
 
 print "done."
 
-ax.arrow(limx_min,limy_min, 10*np.cos(driveangle*np.pi/180),
-		10*np.sin(driveangle*np.pi/180), lw=5.0, color='g')
+ax.arrow(limx_min,limy_min, 10*np.cos(acangle*np.pi/180),
+		10*np.sin(acangle*np.pi/180), lw=2.5, color='blue',
+		head_width=0.15, shape='full')
+ax.arrow(limx_min,limy_min, 10*np.cos(dcangle*np.pi/180),
+		10*np.sin(dcangle*np.pi/180), lw=2.5, color='red',
+		head_width=0.15, shape='full')
+
 
 """
 n = 0
-- 
2.11.4.GIT