From 7be41f2ae4b3d2740e86241118d7b0d049c96dcc Mon Sep 17 00:00:00 2001 From: Rainer Wittmaack Date: Mon, 5 Oct 2009 01:30:23 +0200 Subject: [PATCH] * fixed a bug that caused a QMultiMap to be read out backwards * fixed deterministic diffusion * some minor tweaks and fixes for almost all mpl scripts --- bifandvel.py | 6 ++-- calculator.cpp | 50 +++++++++++++++++------------- calculator.h | 2 +- defangle.py | 57 ++++++++++++++++++++++++++++------ detdiff.py | 38 ++++++++++------------- dimdialog.cpp | 6 ++-- dimdialog.ui | 14 ++++++--- map.py | 4 +-- pinv.py | 98 ++++++++++++++++++++++++++++++---------------------------- 9 files changed, 164 insertions(+), 111 deletions(-) diff --git a/bifandvel.py b/bifandvel.py index ca5547b..98c8f38 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 = '/home/raw/rk_bifurcate.dat' -DATAV = '/home/raw/rk_bif_avg_velo.dat' +DATAX = '../rk_bifurcate.dat' +DATAV = '../rk_bif_avg_velo.dat' FOUT = '../test.png' # some macros @@ -31,7 +31,7 @@ ytkl1 = ["0","L/2","L"] ylim1 = [0-ypad,L+ypad] #ylim2 = [-1-ypad,1+ypad] xlbl1 = r'' -xlbl2 = r'$l$' +xlbl2 = r'$\eta_1$' ylbl1 = r'$x(k\tau)$ mod $L$' ylbl2 = r'$v_x$' res = (1440/80,900/80) # default dpi is 80 diff --git a/calculator.cpp b/calculator.cpp index d8a1e1e..7466385 100644 --- a/calculator.cpp +++ b/calculator.cpp @@ -295,6 +295,7 @@ void Calculator::write_to_file(const QString & name, x.next(); y.next(); phi.next(); + //qDebug() << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value(); out << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value() << endl; } } @@ -308,21 +309,23 @@ void Calculator::write_to_file(const QString & name, //----------------------------------------------------------------------------- -Doub Calculator::mean(VecDoub_I &msd, const Int & samples) +void Calculator::diffusion(VecDoub_I &msd, const Int & samples, Doub & sigma) { - Doub mu,sum,sigma; + Doub mu,sum; VecDoub var(samples,0.0); sum = 0; for(Int i=0;ivaluesMap.value(8); xdot0 = c->valuesMap.value(10); @@ -782,17 +784,17 @@ void Calculator::calc2(const Int n) } break; case 4: - xpi[z][k] = ABS(fmod(yout[0],(2*PI))); + xpi[z][k] = yout[0]; //qDebug() << "xpi["<= s*fromt) { - rkValues_x.insert(r,xpi[z][k]); - rkValues_y.insert(r,ypi[z][k]); - rkValues_phi.insert(r,yout[2]); + rkValues_x.insert(r,ABS(fmod(xpi[z][k],L))); + rkValues_y.insert(r,ABS(fmod(ypi[z][k],L))); + rkValues_phi.insert(r,ABS(fmod(phipi[z][k],L))); rkValues_xdot.insert(r,dydx[0]); rkValues_ydot.insert(r,dydx[1]); rkValues_phidot.insert(r,dydx[2]); @@ -817,6 +819,7 @@ void Calculator::calc2(const Int n) break; case 4: m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L; + avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k; // makes only sense with r = params[10] trajectory_length = sqrt(SQR(yout[0]-x0)+SQR(yout[1])-y0); if (abs(yout[0]-x0) > trajectory_length) { @@ -826,11 +829,12 @@ void Calculator::calc2(const Int n) } bias_angle = params[10]; deflection = (180/PI)*ABS(velocity_angle - bias_angle); + /* qDebug() << velocity_angle << "("<=0;k--) { rkVariance_x.insert(r,sigma_x[k]); rkVariance_y.insert(r,sigma_y[k]); rkVariance_phi.insert(r,sigma_phi[k]); diff --git a/calculator.h b/calculator.h index 788fe76..0392c3b 100644 --- a/calculator.h +++ b/calculator.h @@ -53,7 +53,7 @@ private: void euler(VecDoub_I &y, VecDoub_I &dydx, const Doub x, const Doub h, VecDoub_O &yout); void IC_setup(VecDoub &y, Doub x0, Doub y0, Doub xdot0, Doub ydot0, Doub phi0, Doub params[]); - Doub mean(VecDoub_I & msd, const Int & samples); + void diffusion(VecDoub_I & msd, const Int & samples, Doub & sigma); float gaussian_noise(); void gen_gnuplot(const QString & name); void filter(QMap & map, const Doub & eps); diff --git a/defangle.py b/defangle.py index 0ee3c4b..ce98b03 100644 --- a/defangle.py +++ b/defangle.py @@ -8,15 +8,19 @@ import random from pylab import * from matplotlib.pyplot import * -ANGLE = "../rk_defangle.dat" -VELO = "../rk_bif_avg_velo.dat" +ANGLE = "../data/rk_defangle.dat" +VELO = "../data/rk_defangle_velo.dat" FOUT = "../defangle.png" xlabel = r'$h$' ylabel = r'$\delta y$' res = (1440/80,900/80) -pad1 = 0.01 -pad2 = 0.01 +pad1 = 1 +pad2 = 0.1 +tsize = 16 +lsize = 24 +max = 0.875 +acangle = 58 ###################################################################### @@ -38,15 +42,23 @@ def mean_angle(param,angle): if (param[j] == unique_param[i]): list_angles[i].append(angle[j]) #print len(list_angles[i]) - ret_avg.append(float(sum(list_angles[i])) / len(list_angles[i])) + ret_avg.append(abs(float(sum(list_angles[i])) / + len(list_angles[i]))) #print len(ret_avg) return ret_avg +left, width = 0.1, 0.8 +c = max/2 +rect1 = [left, c+0.1, width, c] +rect2 = [left, 0.1, width, c] + +rad_ac = acangle*(np.pi/180) + fig = plt.figure(figsize=res) -ax1 = fig.add_subplot(211) -ax2 = fig.add_subplot(212) +ax1 = fig.add_axes(rect1) +ax2 = fig.add_axes(rect2) avg_angles = mean_angle(r,psi) avg_velo = mean_angle(r,vx) @@ -55,10 +67,37 @@ unique_r = unique(r) ax1.plot(unique_r, avg_angles, marker='.', color='black') ax2.plot(unique_r, avg_velo, marker='.', color='black') +xticks = np.arange(0,360+1,45) +yticks1 = np.arange(135,315+1,45) +yticks2 = np.arange(0,1.1,0.5) + ax1.set_xlim(r.min()-pad1,r.max()+pad1) -ax1.set_ylim(psi.min()-pad1,psi.max()+pad1) +ax1.set_ylim(yticks1.min(),yticks1.max()) ax2.set_xlim(r.min()-pad2,r.max()+pad2) -ax2.set_ylim(vx.min()-pad2,vx.max()+pad2) +ax2.set_ylim(0,vx.max()+pad2) +ax1.set_xticks([]) +ax2.set_xticks(xticks) +ax1.set_yticks(yticks1) +ax2.set_yticks(yticks2) +ax1.axvline(acangle, color='black', ls=":") +ax1.axvline(acangle-45, color='black', ls="-.") +ax1.axvline(acangle+45, color='black', ls="-.") +ax1.axvline(acangle+180, color='black', ls=":") +ax1.axvline(acangle+180-45, color='black', ls="-.") +ax1.axvline(acangle+180+45, color='black', ls="-.") +ax2.axvline(acangle, color='black', ls=":") +ax2.axvline(acangle+180, color='black', ls=":") +ax2.axvline(acangle+180-45, color='black', ls="-.") +ax2.axvline(acangle+180+45, color='black', ls="-.") +ax2.axvline(acangle-45, color='black', ls="-.") +ax2.axvline(acangle+45, color='black', ls="-.") + +ax1.set_ylabel(r'$\psi$', size=lsize) +ax2.set_ylabel(r'$v_x$', size=lsize) +ax2.set_xlabel(r'$\phi$', size=lsize) +ax2.set_xticklabels(ax2.get_xticks(), size=tsize) +ax1.set_yticklabels(ax1.get_yticks(), size=tsize) +ax2.set_yticklabels(ax2.get_yticks(), size=tsize) plt.savefig(FOUT) diff --git a/detdiff.py b/detdiff.py index 2182f31..9f6da20 100644 --- a/detdiff.py +++ b/detdiff.py @@ -15,9 +15,9 @@ L = 2*np.pi; tau = L/0.1; samples = 3 -xlabel1 = r'$t$' -ylabel1 = r'$$' -xlabel2 = r'$r$' +xlabel1 = r'$k\,[\tau]$' +ylabel1 = r'$\langle {x(k\tau)}^2 - \langle x(k\tau) \rangle^2\rangle$' +xlabel2 = r'$k\,[\tau]$' ylabel2 = r'$D$' res = (1440/80,900/80) # default dpi is 80 @@ -53,8 +53,6 @@ tra_y = get_trajectory(r,y) print "done." t = range(len(tra_x[0])) -print len(tra_x[0]) -print len(x) print t fig = plt.figure(figsize=res) @@ -74,22 +72,20 @@ D_y =[ [] for DUMMYVAR in range(len(unique(r))) ] for i in range(len(unique(r))): - for j in range(len(tra_x[i])): - # if (j > 0): - tmp = tra_x[i][j]/(2*j*tau) + for j in t: + tmp = tra_x[i][j]/(2*(j+1)*tau) + #print tmp,"=",tra_x[i][j],"/",2*(j+1)*tau D_x[i].append(tmp) - + + for i in range(len(unique(r))): - for j in range(len(tra_y[i])): - # if (j > 0): - tmp = tra_y[i][j]/(2*j*tau) + for j in t: + tmp = tra_y[i][j]/(2*(j+1)*tau) D_y[i].append(tmp) -print len(D_x) -print len(D_y) +#print tra_x +#print D_x -print t -print D_x[i] for i in range(len(unique(r))): ax1.plot(t,tra_x[i]) @@ -98,14 +94,14 @@ for i in range(len(unique(r))): ax2.plot(t,D_x[i]) # ax2.plot(t,D_y[i]) -ax1.set_xlim(0,len(tra_x[0])) -#ax1.set_ylim(x.min(),x.max()) +#ax1.set_xlim(0,len(tra_x[0])) +ax1.set_ylim(0,14001) ax1.set_xlabel(xlabel1) ax1.set_ylabel(ylabel1) -ax2.set_xlim(0,len(tra_x[0])) -#ax2.set_ylim(D.min(),D.max()) +#ax2.set_xlim(0,len(tra_x[0])) +ax2.set_ylim(0,2) ax2.set_xlabel(xlabel2) ax2.set_ylabel(ylabel2) plt.savefig(FOUT) -plt.show() +#plt.show() diff --git a/dimdialog.cpp b/dimdialog.cpp index a6aba4c..54b2875 100644 --- a/dimdialog.cpp +++ b/dimdialog.cpp @@ -249,8 +249,8 @@ void DimDialog::on_comboBox_currentIndexChanged() m_ui->veloxSpinBox->setValue(0); m_ui->veloySpinBox->setValue(0); m_ui->phiSpinBox->setValue(0); - m_ui->eta1SpinBox->setValue(1.0); - m_ui->eta2SpinBox->setValue(1.5); + m_ui->eta1SpinBox->setValue(1.5); + m_ui->eta2SpinBox->setValue(1.0); m_ui->FSpinBox->setValue(0.00); m_ui->wSpinBox->setValue(0.1); m_ui->fromr1SpinBox->setValue(0.0); @@ -261,7 +261,7 @@ void DimDialog::on_comboBox_currentIndexChanged() m_ui->tempSpinBox->setValue(0.0); m_ui->aSpinBox->setValue(1.15); m_ui->bSpinBox->setValue(58.0); - m_ui->lSpinBox->setValue(3.0); + m_ui->lSpinBox->setValue(2.5); m_ui->thetaSpinBox->setValue(0.0); m_ui->psiSpinBox->setValue(0.0); m_ui->param1Edit->setText("A"); diff --git a/dimdialog.ui b/dimdialog.ui index ef2069e..e9cf2b0 100644 --- a/dimdialog.ui +++ b/dimdialog.ui @@ -1078,9 +1078,9 @@ - 160 + 150 110 - 61 + 51 20 @@ -1231,12 +1231,18 @@ - 220 + 210 110 - 51 + 61 21 + + 1000 + + + 10 + diff --git a/map.py b/map.py index 1dcbbb1..988a5cb 100644 --- a/map.py +++ b/map.py @@ -21,9 +21,9 @@ xtfreq = 15.2 ytfreq = 1 xpad = 0.001 ypad = 0.001 -ms = 20 +ms = 400 ystep = 0.1 -xstep = 1 +xstep = 0.1 res = (1440/80,900/80) # default dpi is 80 method = 'direct' # mean, direct diff --git a/pinv.py b/pinv.py index a01d0f2..dd3e757 100644 --- a/pinv.py +++ b/pinv.py @@ -20,7 +20,7 @@ tau = L/w evry = 1 # plot periodicity, don't set < 10 (bug) offx = -98 # offset of (x,y) in L offy = 94 -tfrom = 90 # from t +tfrom = 80 # from t tto = 100 # to t acangle = 58 # angle of ac drive dcangle = 0 # angle of dc bias @@ -36,6 +36,8 @@ ylabel = r'$y$' cblabel = r'$\cos(x)\cos(y)+\cos(x)+\cos(y)$' fs_ticks = 16 fs_labels = 24 +cof_only = True +has_dc = False title = 'F = 0.1, alpha = 83, l = 3.0, a = 1.5, eta1 = 1.0, eta2 = 1.5, theta = 0' ################################################################################ @@ -79,7 +81,6 @@ xmin = np.ceil(min(trunc_cofx)/L) ymax = np.ceil(max(trunc_cofy)/L) ymin = np.ceil(min(trunc_cofy)/L) -print xmin*L, xmax*L xdiff = np.ceil(abs(max(trunc_cofx)-min(trunc_cofx))/L) ydiff = np.ceil(abs(max(trunc_cofy)-min(trunc_cofy))/L) @@ -90,17 +91,17 @@ print xdiff, ydiff print xdiff-(xdiff-sper) if xmax > 0: - limx_max = xmax*L - limx_min = (xmax-sper)*L + limx_min = (xmin+2)*L + limx_max = (xmin+2+sper)*L if xmax < 0: - limx_min = xmin*L - limx_max = (xmin+sper)*L + limx_min = (xmax-1-sper)*L + limx_max = (xmax-1)*L if ymax > 0: - limy_max = ymax*L - limy_min = (ymax-sper)*L + limy_min = (ymin-2)*L + limy_max = (ymin-2+sper)*L if ymax < 0: - limy_min = ymin*L - limy_max = (ymin+sper)*L + limy_min = (ymax-sper)*L + limy_max = (ymax)*L print (limx_max-limx_min)/L @@ -142,46 +143,41 @@ cb = plt.colorbar(pot) cb.set_label(cblabel) print "plotting dimer onto the basemap..." -cof = ax.scatter(trunc_cofx, trunc_cofy, c='yellow', marker='o', lw=0, s=msize/4, label='cof') -x1 = ax.scatter(trunc_p1x, trunc_p1y, c='blue', marker='o', lw=0, s=msize, label='x1') -x2 = ax.scatter(trunc_p2x, trunc_p2y, c='purple', marker='o', lw=0, s=msize, label='x2') -#bond = ax.arrow(trunc_p1x, trunc_p1y, trunc_pxdiff, trunc_pydiff, color='r', lw=0.1, label='bond') - - -k=0 -for i in range(len(trunc_t)): - #bond = ax.arrow(p1x[i], p1y[i], pxdiff[i], pydiff[i], color='r', lw=0.1, label='bond') - if trunc_t[i] >= k*every: - #print "plotting because",t[i],"=",k,"*",every,"(",k*every,")" - bond = ax.arrow(trunc_p1x[i], trunc_p1y[i], - trunc_pxdiff[i], trunc_pydiff[i], - color='r', lw=0.1, label='bond', - alpha=0.5) - k = k+1 - elif trunc_t[i] < k*every: - #print "not plotting because",t[i],"!=",k,"*",every,"(",k*every,")" - continue +if (cof_only is False): + cof = ax.scatter(trunc_cofx, trunc_cofy, c='yellow', marker='o', lw=0, s=msize/4, label='cof') + x1 = ax.scatter(trunc_p1x, trunc_p1y, c='blue', marker='o', lw=0, s=msize, label='x1') + x2 = ax.scatter(trunc_p2x, trunc_p2y, c='purple', marker='o', lw=0, s=msize, label='x2') + + k=0 + for i in range(len(trunc_t)): + if trunc_t[i] >= k*every: + #print "plotting because",t[i],"=",k,"*",every,"(",k*every,")" + bond = ax.arrow(trunc_p1x[i], trunc_p1y[i], + trunc_pxdiff[i], + trunc_pydiff[i], color='r', + lw=0.1, label='bond', + alpha=0.5) + k = k+1 + elif trunc_t[i] < k*every: + #print "not plotting because",t[i],"!=",k,"*",every,"(",k*every,")" + continue +elif (cof_only is True): + cof = ax.scatter(trunc_cofx, trunc_cofy, c='red', marker='o', + lw=0, s=msize/4, label='cof') +else: + print "fatal error!" + print "done." 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 -for i in range(len(x)): - if x[i] >= n*L: - s1 = ax.axhline(y=n/2) - n = n+1 - elif x[i] < n*L: - continue -""" +if (has_dc is True): + 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') ax.set_xlabel(xlabel, size=fs_labels) ax.set_ylabel(ylabel, size=fs_labels) @@ -192,10 +188,18 @@ ax.set_ylim(limy_min,limy_max) ax.set_xticklabels(tlabel, size=fs_ticks) ax.set_yticklabels(tlabel, size=fs_ticks) -legend = fig.legend((cof, x1, x2, bond), (r'$\vec{r}_s$', r'$\vec{r}_1$', - r'$\vec{r}_2$', r'$|\vec{r}_1-\vec{r}_2|$'), shadow=True, - fancybox=True) -legend.get_frame().set_alpha(0.75) +if (cof_only is False): + legend = fig.legend((cof, x1, x2, bond), (r'$\vec{r}_s$', r'$\vec{r}_1$', + r'$\vec{r}_2$', r'$|\vec{r}_1-\vec{r}_2|$'), shadow=True, + fancybox=True) +""" +elif (cof_only is True): + legend = fig.legend((cof),(r'$\vec{r}_s$'), shadow=True, fancybox=True) +else: + print "fatal error!" + +#legend.get_frame().set_alpha(0.75) +""" #plt.title(title) -- 2.11.4.GIT