summaryloggraphiclog1graphiclog2commitcommitdifftreerefseditfork
rawpatchinlineside by side (parent: 13a3165)
rawpatchinlineside by side (parent: 13a3165)
| author | Rainer Wittmaack <[email protected]> | |
| Sun, 27 Sep 2009 22:49:25 +0000 (28 00:49 +0200) | ||
| committer | Rainer Wittmaack <[email protected]> | |
| Sun, 27 Sep 2009 22:49:25 +0000 (28 00:49 +0200) |
* you can now rotate U
* converge.py -- plots converge mode data
* defangle.py -- plots deflection angle
* added deflection angle (trajectory angle minus dc bias angle)
* various bugfixes and enhancements
* converge.py -- plots converge mode data
* defangle.py -- plots deflection angle
* added deflection angle (trajectory angle minus dc bias angle)
* various bugfixes and enhancements
diff --git a/bifandvel.py b/bifandvel.py
--- a/bifandvel.py
+++ b/bifandvel.py
import matplotlib.axes as axe
import matplotlib.text as txt
-DATAX = 'rk_bifurcate.dat'
-DATAV = 'rk_bif_avg_velo.dat'
+DATAX = '../rk_bifurcate.dat'
+DATAV = '../rk_bif_avg_velo.dat'
+FOUT = '../test.png'
# some macros
PI = math.pi
# options
xpad = 0.002
-ypad = 0.2
-max = 0.84
-tsize = 12
+ypad = 0.075
+max = 0.875
+tsize = 24
+ticksize = 14
# default values
xtks1 = []
ytks2 = [-1,-0.5,0,0.5,1]
xtkl1 = []
ytkl1 = ["0","L/2","L"]
-ylim1 = [0-0.2,L+0.2]
+ylim1 = [0-ypad,L+ypad]
ylim2 = [-1-ypad,1+ypad]
xlbl1 = r''
xlbl2 = r'$l$'
# draw figure
fig = plt.figure(figsize=res)
-fig.suptitle(title, size=tsize)
+#fig.suptitle(title, size=tsize)
# draw first subplot and setup axis
ax1 = fig.add_axes(rect1)
ax1.set_ylim(ylim1)
ax1.set_yticks(ytks1)
ax1.set_ylabel(ylbl1, size=tsize)
-ax1.set_xticklabels(xtkl1, size=tsize)
-ax1.set_yticklabels(ytkl1, size=tsize)
+ax1.set_xticklabels(xtkl1, size=ticksize)
+ax1.set_yticklabels(ytkl1, size=ticksize)
# put <v> into the plot
ax2.plot(s.r, s.x, c='b', marker=',', ms=0.1, ls='')
ax2.set_yticks(ytks2)
ax2.set_ylabel(ylbl2, size=tsize)
ax2.set_xlabel(xlbl2, size=tsize)
-ax2.set_xticklabels(ax2.get_xticks(), size=tsize)
-ax2.set_yticklabels(ax2.get_yticks(), size=tsize)
+ax2.set_xticklabels(ax2.get_xticks(), size=ticksize)
+ax2.set_yticklabels(ax2.get_yticks(), size=ticksize)
# adjust subplots
plt.subplots_adjust(hspace=0)
# stream the whole mess into a file
-plt.savefig('test.png')
+plt.savefig(FOUT)
-# EOF
+# EOf
diff --git a/calculator.cpp b/calculator.cpp
--- a/calculator.cpp
+++ b/calculator.cpp
file_rk_variance = "../rk_variance.dat";
file_rk_distance = "../rk_msd.dat";
+ file_rk_defangle = "../rk_defangle.dat";
+
+ file_rk_dualbif = "../rk_dualbif.dat";
+
v = c->valuesMap.value(1);
}
//params[14] = e3;
//params[15] = fx;
//params[16] = fy;
+ //params[17] = psi;
+ //params[18] = U;
//static Doub G, xi, xi1, xi2, c, s, st, ax, ay, e0, e1, e2, e3, eta, l1, l2;
- static Doub xi, G, s, c, st, ax, ay;
- static Doub x1, y1, x2, y2, F1x, F1y, F2x, F2y, fx, fy, stx, sty;
+ static Doub xi, G, s, c, st, ax, ay, cp, sp, m1, m2, p1, p2;
+ static Doub x1, y1, x2, y2, F1x, F1y, F2x, F2y, stx, sty, xix1, xix2, xiy1, xiy2;
//for (Int i = 0;i<10;i++)
// qDebug() << "params["<<i<<"] :" << params[i];
break;
case 4:
- //xi1 = gauss.dev();
- //xi2 = gauss.dev();
+ xix1 = gauss.dev();
+ xix2 = gauss.dev();
+ xiy1 = gauss.dev();
+ xiy2 = gauss.dev();
/*
eta = params[5]+params[6];
//qDebug() << y[0] - params[8]*c << "=" << y[0] << "-" << params[8] << "*" << c;
//qDebug() << y[1] - params[8]*c << "=" << y[1] << "-" << params[8] << "*" << s;
- F1x = sin(x1)*cos(y1);//+1.6*sin(x1)+1.0*cos(y1);
- F1y = cos(x1)*sin(y1);//+1.6*cos(x1)+1.0*sin(y1);
- F2x = sin(x2)*cos(y2);//+1.6*sin(x2)+1.0*cos(y1);
- F2y = cos(x2)*sin(y2);//+1.6*cos(x2)+1.0*sin(y2);
+ cp = cos(params[17]);
+ sp = sin(params[17]);
+
+ m1 = cp*x1-sp*y1;
+ p1 = sp*x1+cp*y1;
+ m2 = cp*x2-sp*y2;
+ p2 = sp*x2+cp*y2;
+
+ F1x = params[18]*(sin(m1)*cos(p1)*cp+cos(m1)*sin(p1)*sp+sin(m1)*cp+sin(p1)*sp);//+1.6*sin(x1)+1.0*cos(y1);
+ F1y = params[18]*(-sin(m1)*cos(p1)*sp+cos(m1)*sin(p1)*cp-sin(m1)*sp+sin(p1)*cp);//+1.6*cos(x1)+1.0*sin(y1);
+ F2x = params[18]*(sin(m2)*cos(p2)*cp+cos(m2)*sin(p2)*sp+sin(m2)*cp+sin(p2)*sp);//+1.6*sin(x2)+1.0*cos(y2);
+ F2y = params[18]*(-sin(m2)*cos(p2)*sp+cos(m2)*sin(p2)*cp-sin(m2)*sp+sin(p2)*cp);//+1.6*cos(x2)+1.0*sin(y2);
- dydx[0] = (F1x+F2x+2.0*stx)*params[11];
- dydx[1] = (F1y+F2y+2.0*sty)*params[11];
- dydx[2] = params[12]*(c*F1y-s*F1x)-params[13]*(c*F2y-s*F2x)+params[14]*(c*sty-s*stx);
+ dydx[0] = (F1x+F2x+2.0*stx+sqrt(2*params[9])*(xix1+xix2))*params[11];
+ dydx[1] = (F1y+F2y+2.0*sty+sqrt(2*params[9])*(xiy1+xiy2))*params[11];
+ dydx[2] = params[12]*(c*F1y-s*F1x)-params[13]*(c*F2y-s*F2x)+params[14]*(c*sty-s*stx)+sqrt(2*params[9])*(params[12]*(-s*xix1+c*xiy1)-params[13]*(-s*xix2+c*xiy2));
/*
c = cos(y[2]);
}
break;
case 4:
- while (x.hasNext() && y.hasNext() && phi.hasNext()) {
- x.next();
- y.next();
- phi.next();
- out << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value() << endl;
+ if (map_x.size() > 0 && map_y.size() == 0 && map_phi.size() == 0) {
+ while (x.hasNext()) {
+ x.next();
+ out << x.key() << "\t" << x.value() << endl;
+ }
+ }
+ else if (map_x.size() > 0 && map_y.size() > 0 && map_phi.size() > 0) {
+ while (x.hasNext() && y.hasNext() && phi.hasNext()) {
+ x.next();
+ y.next();
+ phi.next();
+ out << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value() << endl;
+ }
}
break;
}
Doub mu,sum,sigma;
VecDoub var(samples,0.0);
- assert(msd.size() == samples);
sum = 0;
for(Int i=0;i<samples;i++) {
sum += msd[i];
for(Int i=0;i<samples;i++) {
var[i] = SQR(msd[i])-SQR(mu);
}
- assert(var.size() == samples) ;
+
sum = 0;
for(Int j=0;j<samples;j++) {
sum += var[j];
}
}
-void Calculator::IC_setup(VecDoub &y, Doub x0, Doub xdot0, Doub y0, Doub ydot0,
+void Calculator::IC_setup(VecDoub &y, Doub x0, Doub y0, Doub xdot0, Doub ydot0,
Doub phi0, Doub params[])
{
Int n = y.size();
params[6] = c->valuesMap.value(26); // eta2;
params[9] = c->valuesMap.value(32); // T;
params[10] = fmod(c->valuesMap.value(40),360)*PI/180; // theta;
+ params[17] = fmod(c->valuesMap.value(41),360)*PI/180; // psi;
+ params[18] = c->valuesMap.value(42); // U;
Doub e0 = 1.0/(params[5]+params[6]);
Doub e1 = 1.0/(params[5]*params[3]);
Doub e2 = 1.0/(params[6]*params[3]);
{
// setting up a few start values
qDebug() << "calc2 has been called: " << "calc2("<<n<<")";
- Int per;
Normaldev gauss(0,1,time(NULL));
Ran ran(time(NULL));
Doub x1,x2,y1,y2,xtest;
QString param2 = c->str2;
qDebug() << "System" << v;
- Doub x;
+ Doub x, hnew;
+ Doub trajectory_length, deflection;
+ Int k;
Doub xpi[samples][s];
Doub ypi[samples][s];
qDebug() << "number of steps:" << s;
x = t0;
IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
- per = 0;
+ k = 0;
L = 2*PI;
tau = 2*PI/params[4];
- while(per < s) {
+ while(k < s) {
xtest = x+h;
- if (xtest >= per*tau) {
- per++;
+ if (xtest >= k*tau) {
+ k++;
}
derivs(x,y,dydx,params,gauss); // sets first dydx[] for rk4()
rk4(y,dydx,x,h,yout,params,gauss);
- if (per >= s*fromt) {
+ if (k >= s*fromt) {
switch(v) {
case 1:
rkValues_x.insert(x,y[0]);
rkValues_y.insert(x,y[1]);
//rkValues_x.insert(x,xpi);
//rkValues_y.insert(x,ypi);
- rkValues_phi.insert(x,y[2]);
+ rkValues_phi.insert(x,(180/PI)*ABS(fmod(y[2],(2*PI))));
x1 = y[0] + params[7]*cos(y[2]);
y1 = y[1] + params[7]*sin(y[2]);
x2 = y[0] - params[8]*cos(y[2]);
qDebug() << "converge mode";
Int iter;
Doub hh, ss;
- IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
- for (iter=1;(iter*from)<to;iter++) {
- hh = (iter*from);
- ss = (s/hh);
- x = t0;
- qDebug() << "stepsize is" << hh << "with" << ss << "total steps.";
- IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
- for (Int i=0;i<ss&&x<(ss-0.5);i++) {
- derivs(x,y,dydx,params,gauss);
- rk4(y,dydx,x,hh,yout,params,gauss);
- for (Int j=0;j<n;j++)
- y[j] = yout[j];
- x += hh;
- }
- switch(v) {
- case 1:
- rkValues_x.insert(hh,yout[0]);
- break;
- case 4:
- rkValues_x.insert(hh,yout[0]);
- rkValues_y.insert(hh,yout[1]);
- rkValues_phi.insert(hh,yout[2]);
- break;
- }
- x = t0;
- IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
- for (Int i=0;i<ss&&x<(ss-0.5);i++) {
- derivs(x,y,dydx,params,gauss);
- euler(y,dydx,x,hh,yout);
- for (Int j=0;j<n;j++)
- y[j] = yout[j];
- x += hh;
- }
- switch(v) {
- case 1:
- euValues_x.insert(hh,yout[0]);
- break;
- case 4:
- euValues_x.insert(hh,yout[0]);
- euValues_y.insert(hh,yout[1]);
- euValues_phi.insert(hh,yout[2]);
- break;
- }
- }
+ for (iter=1;(iter*from)<=to;iter++) {
+ for (Int z=0;z<samples;z++) {
+ qDebug() << "at step" << z << "of" << samples;
+ x = t0;
+ hh = (iter*from);
+ ss = (s/hh);
+ qDebug() << "stepsize is" << hh << "with" << ss << "total steps.";
+ //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
+ IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
+ k = 0;
+ L = 2*PI;
+ tau = 2*PI/params[4];
+ while(k < s) {
+ //qDebug() << k << "<" << s;
+ //qDebug() << k << tau << k*tau;
+ xtest = x+hh;
+ if (xtest >= k*tau) {
+ //qDebug() << xtest << ">=" << k << "*" << tau << "("<<k*tau<<")";
+ hnew = k*tau-x;
+ derivs(x,y,dydx,params,gauss);
+ rk4(y,dydx,x,hnew,yout,params,gauss);
+ k++;
+ x += hnew;
+ //qDebug() << x;
+ } else {
+ //qDebug() << xtest << "=" << x<<"+"<<hh;
+ derivs(x,y,dydx,params,gauss);
+ rk4(y,dydx,x,hh,yout,params,gauss); // xold
+ x += hh;
+ }
+ for (Int j=0;j<n;j++)
+ y[j] = yout[j];
+ } // end while
+ switch(v) {
+ case 1:
+ rkValues_x.insert(hh,yout[0]);
+ break;
+ case 4:
+ rkValues_x.insert(hh,yout[0]);
+ rkValues_y.insert(hh,yout[1]);
+ rkValues_phi.insert(hh,yout[2]);
+ break;
+ }
+ x = t0;
+ //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
+ IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
+ k = 0;
+ L = 2*PI;
+ tau = 2*PI/params[4];
+ while (k < s) {
+ xtest = x+hh;
+ if (xtest >= k*tau) {
+ hnew = k*tau-x;
+ derivs(x,y,dydx,params,gauss);
+ euler(y,dydx,x,hnew,yout);
+ k++;
+ x += hnew;
+ } else {
+ derivs(x,y,dydx,params,gauss);
+ euler(y,dydx,x,hh,yout);
+ x += hh;
+ }
+ for (Int j=0;j<n;j++)
+ y[j] = yout[j];
+ } // end while
+ switch(v) {
+ case 1:
+ euValues_x.insert(hh,yout[0]);
+ break;
+ case 4:
+ euValues_x.insert(hh,yout[0]);
+ euValues_y.insert(hh,yout[1]);
+ euValues_phi.insert(hh,yout[2]);
+ break;
+ } // end switch
+ } // end samples
+ } // end for
qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
write_to_file(file_rk_conv, rkValues_x, rkValues_y, rkValues_phi, params);
write_to_file(file_eu_conv, euValues_x, euValues_y, euValues_phi, params);
//for(Int i = 0;i<=10;i++)
// qDebug() << params[i];
// inner for: iterate over runge-kutta
- Int k = 0;
+ k = 0;
Doub hnew;
L = 2*PI;
tau = 2*PI/params[4];
break;
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)));
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);
break;
} // end switch
} // samples loop
VecDoub sigma_phi(s,0.0);
//qDebug() << "samples:" << samples << "steps:" << s;
for(Int i=0;i<s;i++) {
- //qDebug() << "["<<i<<"]:-----------------------------------------";
+ //qDebug() << "["<<i<<"]:---------------------------------";
for(Int j=0;j<samples;j++) {
//qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
msd_x[j] = xpi[j][i];
write_to_file(file_rk_avg, rkAvgvelo_x, rkAvgvelo_y, rkAvgvelo_phi, params);
write_to_file(file_rk_variance, rkVariance_x, rkVariance_y, rkVariance_phi,params);
write_to_file(file_rk_distance, rkDist_x, rkDist_y, rkDist_phi, params);
+ write_to_file(file_rk_defangle, rkDefAngle, NullMap, NullMap, params);
break;
case 3:
qDebug() << "dual bifurcation mode";
qDebug() << "Number of samples : " << samples;
Doub stp1, stp2;
QString vx_color, vy_color;
- QString name = "dual_bif.dat";
+ QString name = "tmp.dat";
QFile file(name);
QTextStream out(&file);
file.open(QIODevice::WriteOnly);
} // end 2. outer for
} // end 1. outer for
file.close();
+ qDebug() << "moving" << name << "to safety";
+ file.copy(name, file_rk_dualbif);
break;
} // end switch
} // end if
if (eul) {
+ Doub eu_h = h/100;
switch(M) {
case 0:
qDebug() << "(euler) standard mode";
- qDebug() << "stepsize:" << h;
+ qDebug() << "stepsize:" << eu_h;
qDebug() << "number of steps:" << s;
x = t0;
IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
- per = 0;
+ k = 0;
L = 2*PI;
- tau = 2*PI/params[4];
- while(per < s) {
- xtest = x+h;
- if (xtest >= per*tau) {
- per++;
- }
- derivs(x,y,dydx,params,gauss); // sets first dydx[] for rk4()
- euler(y,dydx,x,h,yout);
- if (per >= s*fromt) {
- switch(v) {
+ Doub hnew;
+ Doub xnew;
+ tau = 0.1;
+ while(k < s) {
+ xtest = x+eu_h;
+ if (xtest >= k*tau) {
+ hnew = k*tau-x;
+ xnew = x+hnew;
+ derivs(x,y,dydx,params,gauss); // sets first dydx[] for rk4()
+ euler(y,dydx,x,hnew,yout);
+ if (k >= s*fromt) {
+ switch(v) {
case 1:
- euValues_x.insert(x,y[0]);
- break;
+ euValues_x.insert(xnew,yout[0]);
+ break;
case 4:
- //xpi = ABS(fmod(yout[0],(2*PI)));
- //ypi = ABS(fmod(yout[1],(2*PI)));
- euValues_x.insert(x,y[0]);
- euValues_y.insert(x,y[1]);
- //rkValues_x.insert(x,xpi);
- //rkValues_y.insert(x,ypi);
- euValues_phi.insert(x,y[2]);
- x1 = y[0] + params[7]*cos(y[2]);
- y1 = y[1] + params[7]*sin(y[2]);
- x2 = y[0] - params[8]*cos(y[2]);
- y2 = y[1] - params[8]*sin(y[2]);
- euValues_p1x.insert(x,x1);
- euValues_p1y.insert(x,y1);
- euValues_p2x.insert(x,x2);
- euValues_p2y.insert(x,y2);
- break;
+ //xpi = ABS(fmod(yout[0],(2*PI)));
+ //ypi = ABS(fmod(yout[1],(2*PI)));
+ euValues_x.insert(xnew,yout[0]);
+ euValues_y.insert(xnew,yout[1]);
+ //rkValues_x.insert(x,xpi);
+ //rkValues_y.insert(x,ypi);
+ euValues_phi.insert(xnew,yout[2]);
+ x1 = yout[0] + params[7]*cos(yout[2]);
+ y1 = yout[1] + params[7]*sin(yout[2]);
+ x2 = yout[0] - params[8]*cos(yout[2]);
+ y2 = yout[1] - params[8]*sin(yout[2]);
+ euValues_p1x.insert(xnew,x1);
+ euValues_p1y.insert(xnew,y1);
+ euValues_p2x.insert(xnew,x2);
+ euValues_p2y.insert(xnew,y2);
+ break;
+ }
}
+ k++;
}
+ derivs(x,y,dydx,params,gauss);
+ euler(y,dydx,x,eu_h,yout); // xold
for (Int j=0;j<n;j++)
y[j] = yout[j];
- x += h;
+ x += eu_h;
} // end while
qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
tau = 2*PI/params[4];
//qDebug() << "tau: " << tau << "\tL: " << L;
while (k<s) {
- xtest = x+h; // test next x increment
+ xtest = x+eu_h; // test next x increment
//if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
if (xtest >= k*tau) {
//qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
- //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
+ //qDebug() << "new eu_h is" << k << "*" << tau << "-" << x << "=" << hnew;
hnew = k*tau-x;
//qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
derivs(x,y,dydx,params,gauss);
k++;
}
derivs(x,y,dydx,params,gauss);
- euler(y,dydx,x,h,yout); // xold
+ euler(y,dydx,x,eu_h,yout); // xold
for (Int j=0;j<n;j++)
- y[j] = yout[j]; // y[x+h] !!!!!!!
- x += h; // do it for real
+ y[j] = yout[j]; // y[x+eu_h] !!!!!!!
+ x += eu_h; // do it for real
} // end while
//qDebug() << "at" << x << "in time for k = " << k;
switch(v) {
/*
* we need to do something about this
*/
- filter(euValues_x,eps);
- filter(euValues_y,eps);
- filter(euValues_xdot,eps);
- filter(euValues_ydot,eps);
- filter(euAvgvelo_x, eps);
- filter(euAvgvelo_y, eps);
- filter(euAvgvelo_phi,eps);
qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
qDebug() << "sizeof(euValues_xdot) =" << euValues_xdot.size();
tau = 2*PI/params[4];
//qDebug() << "tau: " << tau << "\tL: " << L;
while (k<s) {
- xtest = x+h; // test next x increment
+ xtest = x+eu_h; // test next x increment
//if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
if (xtest >= k*tau) {
//qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
- //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
+ //qDebug() << "new eu_h is" << k << "*" << tau << "-" << x << "=" << hnew;
hnew = k*tau-x;
//qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
derivs(x,y,dydx,params,gauss);
k++;
}
derivs(x,y,dydx,params,gauss);
- euler(y,dydx,x,h,yout); // xold
+ euler(y,dydx,x,eu_h,yout); // xold
for (Int j=0;j<n;j++)
- y[j] = yout[j]; // y[x+h] !!!!!!!
- x += h; // do it for real
+ y[j] = yout[j]; // y[x+eu_h] !!!!!!!
+ x += eu_h; // do it for real
} // end while
Int m_x, m_y, m_phi;
m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
diff --git a/calculator.h b/calculator.h
--- a/calculator.h
+++ b/calculator.h
void rk4(VecDoub_I &y, VecDoub_I &dydx, const Doub x, const Doub h, VecDoub_O &yout,
const Doub params[], Normaldev &gauss);
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 xdot0, Doub y0, Doub ydot0, Doub phi0,
+ 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);
float gaussian_noise();
rkValues_p1y, rkValues_p2y,
rkValues_phi, rkValues_phidot,
rkVariance_x, rkVariance_y, rkVariance_phi,
- rkDist_x, rkDist_y, rkDist_phi;
+ rkDist_x, rkDist_y, rkDist_phi,
+ rkDefAngle;
QMultiMap<double, double> euValues_x, euValues_y,
euValues_xdot, euValues_ydot,
euValues_xtemp, euValues_ytemp, euValues_phitemp,
euValues_p1x, euValues_p2x,
euValues_p1y, euValues_p2y,
- euValues_phi, euValues_phidot;
+ euValues_phi, euValues_phidot,
+ euDefAngle;
+
+ QMultiMap<double, double> NullMap; // dummy var for use with write_to_file()
QString file_rk_bif, file_rk_avg, file_rk_out, file_rk_pt1, file_rk_pt2,
- file_rk_biv, file_rk_conv, file_rk_variance, file_rk_distance;
+ file_rk_biv, file_rk_conv, file_rk_variance, file_rk_distance,
+ file_rk_defangle, file_rk_dualbif;
QString file_eu_bif, file_eu_avg, file_eu_out, file_eu_pt1, file_eu_pt2,
- file_eu_biv, file_eu_conv;
+ file_eu_biv, file_eu_conv, file_eu_defangle, file_eu_dualbif;
Int v;
Doub T;
diff --git a/converge.py b/converge.py
--- /dev/null
+++ b/converge.py
@@ -0,0 +1,94 @@
+import math
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.mlab as mlab
+from matplotlib.pyplot import *
+from pylab import *
+
+DATA_RK = "../rk_conv.dat"
+DATA_EU = "../eu_conv.dat"
+FOUT = "../converge.png"
+
+xlabel = r'$h$'
+ylabel = r'$\delta y$'
+res = (1440/80,900/80) # default dpi is 80
+
+######################################################################
+
+def mean_var(param,var):
+ unique_param = unique(param)
+ list_var = [ [] for DUMMYVAR in range(len(unique_param)) ]
+ ret_avg = []
+ for i in range(len(unique_param)):
+ print "now:", unique_param[i]
+ for j in range(len(var)):
+ if (param[j] == unique_param[i]):
+ list_var[i].append(var[j])
+ #assert len(list_var[i]) == len(param/unique_param)
+ ret_avg.append(float(sum(list_var[i])) /
+ len(list_var[i]))
+ #print list_var[i], ret_avg[i]
+ #assert list_var[i] == ret_avg[i]
+ return ret_avg
+
+print "reading data..."
+
+rk = mlab.csv2rec(DATA_RK, delimiter='\t', comments='#')
+eu = mlab.csv2rec(DATA_EU, delimiter='\t', comments='#')
+
+h_rk = rk.r
+h_eu = eu.r
+x_rk = rk.x
+x_eu = eu.x
+
+h_min = h_rk.min()
+h_max = h_rk.max()
+
+print "calculating mean..."
+
+rk_mean = mean_var(h_rk, x_rk)
+eu_mean = mean_var(h_eu, x_eu)
+
+#print rk_mean
+
+
+fig = plt.figure(figsize=res)
+
+#ax1 = fig.add_subplot(211)
+ax2 = fig.add_subplot(111)
+
+#ax1.plot(unique(h_rk), rk_mean, marker=',', ls='')
+#ax1.plot(unique(h_eu), eu_mean, marker=',', ls='')
+
+#ax1.axhline(x_rk[0], color='r')
+#ax1.axhline(x_eu[0], color='r')
+
+print "calculating error..."
+
+err_rk = []
+err_eu = []
+
+for i in range(len(rk_mean)):
+ tmp = np.abs(rk_mean[0] - rk_mean[i])
+ err_rk.append(tmp)
+
+for i in range(len(eu_mean)):
+ tmp = np.abs(rk_mean[0] - eu_mean[i])
+ err_eu.append(tmp)
+
+print "plotting..."
+3
+ax2.loglog(unique(h_rk), err_rk, marker=',', ls='', ms=0.1)
+ax2.loglog(unique(h_eu), err_eu, marker=',', ls='', ms=0.1)
+
+pad = h_max
+
+ax2.set_xlabel(xlabel)
+ax2.set_ylabel(ylabel)
+ax2.set_xlim(h_min,h_max)
+ax2.set_ylim(min(err_eu),max(err_eu))
+
+#ax1.set_xlim()
+#ax1.set_ylim(-3.3,-3.2)
+plt.savefig(FOUT)
diff --git a/defangle.py b/defangle.py
--- /dev/null
+++ b/defangle.py
@@ -0,0 +1,59 @@
+import matplotlib as mpl
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+import matplotlib.axes as axe
+import matplotlib.collections as collections
+import numpy as np
+import random
+from pylab import *
+from matplotlib.pyplot import *
+
+ANGLE = "../rk_defangle.dat"
+VELO = "../rk_bif_avg_velo.dat"
+FOUT = "../defangle.png"
+
+xlabel = r'$h$'
+ylabel = r'$\delta y$'
+
+######################################################################
+
+angle = mlab.csv2rec(ANGLE, delimiter='\t', comments='#')
+velo = mlab.csv2rec(VELO, delimiter='\t', comments='#')
+
+r = angle.r
+psi = angle.x
+vx = velo.x
+
+
+
+def mean_angle(param,angle):
+ unique_param = unique(param)
+ list_angles = [ [] for DUMMYVAR in range(len(unique_param)) ]
+ ret_avg = []
+ for i in range(len(unique_param)):
+ for j in range(len(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]))
+ #print len(ret_avg)
+ return ret_avg
+
+
+fig = plt.figure()
+
+ax1 = fig.add_subplot(211)
+ax2 = fig.add_subplot(212)
+
+avg_angles = mean_angle(r,psi)
+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.set_xlim(r.min(),r.max())
+ax1.set_ylim(psi.min(),psi.max())
+
+#plt.savefig(FOUT)
+
+plt.show()
diff --git a/detdiff.py b/detdiff.py
--- a/detdiff.py
+++ b/detdiff.py
from pylab import *\r
from matplotlib.pyplot import *\r
\r
-FIN = 'quplot/rk_variance.dat'\r
-FOUT = 'diffusion.png'\r
+FIN = '../rk_variance.dat'\r
+FOUT = '../diffusion.png'\r
\r
L = 2*np.pi;\r
+tau = L/0.1;\r
samples = 3\r
\r
xlabel1 = r'$t$'\r
#return (ret_x, ret_y)\r
return (ret_x)\r
\r
-data = mlab.csv2rec(FIN, delimiter='\t')\r
+print "reading data..."\r
\r
+data = mlab.csv2rec(FIN, comments='#', delimiter='\t')\r
r = data.r\r
x = data.x\r
y = data.y\r
\r
+print "done."\r
+\r
+print "extracting trajectories..."\r
tra_x = get_trajectory(r,x)\r
tra_y = get_trajectory(r,y)\r
+print "done."\r
\r
t = range(len(tra_x[0]))\r
print len(tra_x[0])\r
ax2 = fig.add_subplot(122)\r
\r
\r
-#D =[ [] for DUMMYVAR in range(len(unique(r))) ]\r
-D_x = []\r
-D_y = []\r
+D_x =[ [] for DUMMYVAR in range(len(unique(r))) ]\r
+D_y =[ [] for DUMMYVAR in range(len(unique(r))) ]\r
+\r
+#D_x = []\r
+#D_y = []\r
\r
#dmin = min(unique(D))\r
#dmax = max(unique(D))\r
for i in range(len(unique(r))):\r
for j in range(len(tra_x[i])):\r
# if (j > 0):\r
- tmp = (tra_x[i][j] - tra_x[i][j-1])/2\r
- D_x.append(tmp)\r
+ tmp = tra_x[i][j]/(2*j*tau)\r
+ D_x[i].append(tmp)\r
\r
for i in range(len(unique(r))):\r
for j in range(len(tra_y[i])):\r
# if (j > 0):\r
- tmp = (tra_y[i][j] - tra_y[i][j-1])/2\r
- D_y.append(tmp)\r
+ tmp = tra_y[i][j]/(2*j*tau)\r
+ D_y[i].append(tmp)\r
\r
print len(D_x)\r
print len(D_y)\r
\r
+print t\r
+print D_x[i]\r
+\r
for i in range(len(unique(r))):\r
ax1.plot(t,tra_x[i])\r
\r
-ax2.plot(unique(r),D_x)\r
-ax2.plot(unique(r),D_y)\r
+for i in range(len(unique(r))):\r
+ ax2.plot(t,D_x[i])\r
+# ax2.plot(t,D_y[i])\r
\r
ax1.set_xlim(0,len(tra_x[0]))\r
#ax1.set_ylim(x.min(),x.max())\r
ax1.set_xlabel(xlabel1)\r
ax1.set_ylabel(ylabel1)\r
-ax2.set_xlim(r.min(),r.max())\r
+ax2.set_xlim(0,len(tra_x[0]))\r
#ax2.set_ylim(D.min(),D.max())\r
ax2.set_xlabel(xlabel2)\r
ax2.set_ylabel(ylabel2)\r
diff --git a/dimdialog.cpp b/dimdialog.cpp
--- a/dimdialog.cpp
+++ b/dimdialog.cpp
m_ui->fromtSpinBox->setValue(0);
m_ui->tempSpinBox->setValue(0.0);
m_ui->aSpinBox->setValue(1.15);
- m_ui->bSpinBox->setValue(83.0);
+ m_ui->bSpinBox->setValue(58.0);
m_ui->lSpinBox->setValue(3.0);
m_ui->thetaSpinBox->setValue(0.0);
+ m_ui->psiSpinBox->setValue(0.0);
m_ui->param1Edit->setText("A");
m_ui->param2Edit->setText("l");
m_ui->fromr1SpinBox->setValue(0.0);
m_ui->tor1SpinBox->setValue(90.0);
m_ui->fromr2SpinBox->setValue(1.0);
m_ui->tor2SpinBox->setValue(10.0);
+ m_ui->USpinBox->setValue(1.0);
break;
}
}
d->valuesMap.insert(13,m_ui->stepsizeSpinBox->value());
d->valuesMap.insert(14,m_ui->stepsSpinBox->value());
d->valuesMap.insert(40,m_ui->thetaSpinBox->value());
+ d->valuesMap.insert(41,m_ui->psiSpinBox->value());
+ d->valuesMap.insert(42,m_ui->USpinBox->value());
d->valuesMap.insert(29,0); /* disable euler for now */
diff --git a/dimdialog.ui b/dimdialog.ui
--- a/dimdialog.ui
+++ b/dimdialog.ui
<property name="text">
<string>th</string>
</property>
+ <property name="buddy">
+ <cstring>thetaSpinBox</cstring>
+ </property>
+ </widget>
+ <widget class="QDoubleSpinBox" name="psiSpinBox">
+ <property name="geometry">
+ <rect>
+ <x>210</x>
+ <y>90</y>
+ <width>62</width>
+ <height>22</height>
+ </rect>
+ </property>
+ </widget>
+ <widget class="QLabel" name="psiLabel">
+ <property name="geometry">
+ <rect>
+ <x>190</x>
+ <y>90</y>
+ <width>21</width>
+ <height>16</height>
+ </rect>
+ </property>
+ <property name="text">
+ <string>psi</string>
+ </property>
+ <property name="buddy">
+ <cstring>psiSpinBox</cstring>
+ </property>
+ </widget>
+ <widget class="QDoubleSpinBox" name="USpinBox">
+ <property name="geometry">
+ <rect>
+ <x>30</x>
+ <y>120</y>
+ <width>62</width>
+ <height>23</height>
+ </rect>
+ </property>
+ </widget>
+ <widget class="QLabel" name="ULabel">
+ <property name="geometry">
+ <rect>
+ <x>10</x>
+ <y>120</y>
+ <width>21</width>
+ <height>18</height>
+ </rect>
+ </property>
+ <property name="text">
+ <string>U</string>
+ </property>
+ <property name="buddy">
+ <cstring>USpinBox</cstring>
+ </property>
</widget>
</widget>
<widget class="QComboBox" name="comboBox">
<property name="text">
<string>filter</string>
</property>
+ <property name="buddy">
+ <cstring>filterSpinBox</cstring>
+ </property>
</widget>
</widget>
<widget class="QGroupBox" name="groupBox_4">
--- a/map.py
+++ b/map.py
from pylab import *\r
from matplotlib.pyplot import *\r
\r
-FIN = 'dual_bif.dat'\r
-FOUT = 'map_alphaltheta45f01.png'\r
+FIN = 'tmp.dat'\r
+FOUT = '../map_test.png'\r
\r
-samples = 3\r
+samples = 10\r
ticksize = 14\r
labelsize = 24\r
title = "a = 1.15, w = 0.1, F = 0.00, eta1 = 1.0, eta2 = 1.5, theta = 0"\r
-xlabel = r'$\alpha$'\r
-ylabel = r'$|l|$'\r
-tfreq = 10\r
-ms = 40\r
+xlabel = r'$a$'\r
+ylabel = r'$l$'\r
+xtfreq = 15.2\r
+ytfreq = 1\r
+xpad = 0.001\r
+ypad = 0.001\r
+ms = 20\r
+ystep = 0.1\r
+xstep = 1\r
res = (1440/80,900/80) # default dpi is 80\r
method = 'direct' # mean, direct\r
\r
ret_label = r'other rational $v$'\r
return (ret_label)\r
\r
-def assign_color2(colors,xar,yar):\r
- """\r
- number velocity\r
- 0 everything else\r
- 1 v = 0\r
- 2 v = -1.0\r
- 3 v = -0.5\r
- 4 v = 0.5\r
- 5 v = 1.0\r
- 6 v = -0.5 && v = 0.5\r
- 7 v = -1.0 && v = 1.0\r
- """\r
- #print "---------------------------------------------------------"\r
+def assign_color2(colors,xar,yar,clist):\r
+ print "---------------------------------------------------------"\r
unique_colors = unique(colors)\r
- ret_xcolors = [ [] for DUMMYVAR in range(8) ]\r
- ret_ycolors = [ [] for DUMMYVAR in range(8) ]\r
+ ret_xcolors = [ [] for DUMMYVAR in range(len(clist)) ]\r
+ ret_ycolors = [ [] for DUMMYVAR in range(len(clist)) ]\r
i = 0\r
assert len(xar) == len(yar) == len(colors)\r
xarlen_trunc = len(xar)\r
- while (i < xarlen_trunc):\r
- #print len(xar)\r
- #print xarlen_trunc\r
- tmp = []\r
- for j in range(samples):\r
- #print i,j,i+j\r
+ while i < xarlen_trunc:\r
+ tmp = []\r
+ for j in range(samples):\r
assert i+j <= len(xar)\r
- assert yar[i+j] == yar[i] and xar[i+j] == xar[i]\r
- tmp.append(colors[i+j])\r
+ #assert yar[i+j] == yar[i] and xar[i+j] == xar[i]\r
+ tmp.append(colors[i+j])\r
+# print len(tmp), samples\r
assert len(tmp) == samples\r
- if 'green' in tmp and 'orange' in tmp:\r
- ret_xcolors[6].append(xar[i])\r
- ret_ycolors[6].append(yar[i])\r
- elif 'blue' in tmp and 'red' in tmp:\r
- ret_xcolors[7].append(xar[i])\r
- ret_ycolors[7].append(yar[i])\r
- elif len(unique(tmp)) == 1 and 'grey' in tmp:\r
- ret_xcolors[1].append(xar[i])\r
- ret_ycolors[1].append(yar[i])\r
+ if 'green' in tmp and 'orange' in tmp:\r
+ ret_xcolors[6].append(xar[i])\r
+ ret_ycolors[6].append(yar[i])\r
+ elif 'blue' in tmp and 'red' in tmp:\r
+ ret_xcolors[7].append(xar[i])\r
+ ret_ycolors[7].append(yar[i])\r
+ elif len(unique(tmp)) == 1 and 'grey' in tmp:\r
+ ret_xcolors[1].append(xar[i])\r
+ ret_ycolors[1].append(yar[i])\r
elif len(unique(tmp)) == 1 and 'blue' in tmp:\r
ret_xcolors[2].append(xar[i])\r
ret_ycolors[2].append(yar[i])\r
elif len(unique(tmp)) == 1 and 'red' in tmp:\r
ret_xcolors[5].append(xar[i])\r
ret_ycolors[5].append(yar[i])\r
- else:\r
- ret_xcolors[0].append(xar[i])\r
- ret_ycolors[0].append(yar[i])\r
+ else:\r
+ ret_xcolors[0].append(xar[i])\r
+ ret_ycolors[0].append(yar[i])\r
i = i+samples\r
#print "increment:",i\r
- #print "---------------------------------------------------------"\r
+ for j in range(8):\r
+ print clist[j],":",len(ret_xcolors[j])\r
+ print "---------------------------------------------------------"\r
return (ret_xcolors, ret_ycolors)\r
\r
def assign_color(colors,xar,yar):\r
print "---------------------------------------------------------"\r
return (ret_xcolors, ret_ycolors)\r
\r
+print "reading data..."\r
r = mlab.csv2rec(FIN, delimiter='\t', comments='#')\r
\r
x = r.x\r
vy = r.vy\r
\r
assert len(x) == len(y) == len(vx) == len(vy)\r
+print "done."\r
\r
xmin = np.min(unique(x))\r
xmax = np.max(unique(x))\r
ymax = np.max(unique(y))\r
\r
# assumes that the step distance doesn't change over the course of time\r
-ystep = 0.1\r
-xstep = 1\r
rat = xstep/ystep\r
print "size of xstep:", xstep\r
print "size of ystep:", ystep\r
\r
label_list = [\r
r'other rational $v$', # 0 black\r
- r'$v = 0$', # 1 grey\r
- r'$v| = -1$', # 2 blue\r
- r'$v = -0.5$', # 3 green\r
- r'$v = 0.5$', # 4 orange\r
- r'$v = 1.0$', # 5 red\r
- r'$v = \pm 0.5$', # 6 teal\r
- r'$v = \pm 1.0$' # 7 purple\r
+ r'$v = 0$', # 1 grey\r
+ r'$v = -1$', # 2 blue\r
+ r'$v = -0.5$', # 3 green\r
+ r'$v = 0.5$', # 4 orange\r
+ r'$v = 1.0$', # 5 red\r
+ r'$v = \pm 0.5$', # 6 teal\r
+ r'$v = \pm 1.0$' # 7 purple\r
]\r
\r
xpadmin = xmin-xstep\r
ypadmax = ymax+ystep\r
xlim = (xpadmin, xpadmax) # small padding between plot and\r
ylim = (ypadmin, ypadmax) # axis so they don't overlap\r
-xticks = np.linspace(xmin, xmax, tfreq)\r
-yticks = np.linspace(ymin, ymax+ystep, tfreq)\r
+xticks = np.arange(xmin, xmax+xpad, xtfreq)\r
+yticks = np.arange(ymin, ymax+ypad, ytfreq)\r
\r
\r
fig = plt.figure(figsize=res)\r
-fig.suptitle(title, size=labelsize)\r
+#fig.suptitle(title, size=labelsize)\r
\r
max = 0.84\r
c = max/2\r
\r
if (method == 'mean'):\r
\r
+ print "assigning colors...."\r
xcolors_x,xcolors_y = assign_color(x_color,x,y)\r
ycolors_x,ycolors_y = assign_color(y_color,x,y)\r
-\r
+ print "done."\r
+ \r
+ print "plotting..."\r
for k in range(len(unique(x_color))):\r
x_xtli = xcolors_x[k]\r
x_ytli = xcolors_y[k]\r
labelx = assign_label(unique(x_color)[k])\r
ax1.scatter(x_xtli, x_ytli, marker='s', lw=0, s=ms, alpha=0.5,\r
- c = unique(x_color)[k], label=labelx)\r
+ c = unique(x_color)[k], label=labelx,\r
+ edgecolor='face')\r
\r
for k in range(len(unique(y_color))):\r
y_xtli = ycolors_x[k]\r
y_ytli = ycolors_y[k]\r
labely = assign_label(unique(y_color)[k])\r
ax2.scatter(y_xtli, y_ytli, marker='s', lw=0, s=ms, alpha=0.5,\r
- c = unique(y_color)[k], label=labely)\r
+ c = unique(y_color)[k], label=labely,\r
+ edgecolor='none')\r
\r
+ print "done."\r
+ \r
elif (method == 'direct'):\r
\r
- xcolors_x,xcolors_y = assign_color2(x_color,x,y)\r
- ycolors_x,ycolors_y = assign_color2(y_color,x,y)\r
-\r
- for k in range(len(color_list)):\r
- x_xtli = xcolors_x[k]\r
- x_ytli = xcolors_y[k]\r
- #labelx = assign_label(color_list)[k]\r
- ax1.scatter(x_xtli, x_ytli, marker='s', lw=0, s=ms, alpha=0.5,\r
- c = color_list[k], label=label_list[k])\r
-\r
- for k in range(len(color_list)):\r
- y_xtli = ycolors_x[k]\r
- y_ytli = ycolors_y[k]\r
- #labely = assign_label(color_list[k])\r
- ax2.scatter(y_xtli, y_ytli, marker='s', lw=0, s=ms, alpha=0.5,\r
- c = color_list[k], label=label_list[k])\r
-\r
+ print "assigning colors..."\r
+ xcolors_x,xcolors_y = assign_color2(x_color,x,y,color_list)\r
+ ycolors_x,ycolors_y = assign_color2(y_color,x,y,color_list)\r
+ print "done."\r
+\r
+ #for i in range(8):\r
+ # print i, xcolors_x[i],xcolors_y[i];\r
+ \r
+ print "plotting..."\r
+ for k in range(len(xcolors_x)):\r
+ if len(xcolors_x[k]) > 0:\r
+ x_xtli = xcolors_x[k]\r
+ x_ytli = xcolors_y[k]\r
+ ax1.scatter(x_xtli, x_ytli, marker='s', lw=0,\r
+ s=ms, alpha=1, c=color_list[k],\r
+ label=label_list[k])\r
+\r
+ for k in range(len(ycolors_x)):\r
+ if len(ycolors_x[k]) > 0:\r
+ y_xtli = ycolors_x[k]\r
+ y_ytli = ycolors_y[k]\r
+ ax2.scatter(y_xtli, y_ytli, marker='s', lw=0,\r
+ s=ms, alpha=1, c=color_list[k],\r
+ label=label_list[k])\r
+\r
+ print "done."\r
else:\r
print "fatal error: specify a valid method"\r
\r
-leg1 = ax1.legend(loc=0, shadow=True, fancybox=True, scatterpoints=1,\r
- markerscale=1)\r
-leg2 = ax2.legend(loc=0, shadow=True, fancybox=True, scatterpoints=1,\r
- markerscale=1)\r
+leg1 = ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, shadow=True,\r
+ fancybox=True, scatterpoints=1, markerscale=1,\r
+ borderaxespad=0.)\r
+leg2 = ax2.legend(bbox_to_anchor=(1.05, 1), loc=2, shadow=True,\r
+ fancybox=True, scatterpoints=1, markerscale=1,\r
+ borderaxespad=0.)\r
for t in leg1.get_texts():\r
t.set_fontsize('small')\r
for t in leg2.get_texts():\r
leg2.get_frame().set_alpha(0.75)\r
\r
plt.subplots_adjust(hspace=0)\r
+print "flushing into a file..."\r
plt.savefig(FOUT)\r
-plt.show()\r
+print "done."\r
+#plt.show()\r
\r
\r
\r
--- a/pinv.py
+++ b/pinv.py
+#!/usr/bin/python\r
+# -*- coding: iso-8859-1 -*-\r
+\r
import matplotlib as mpl\r
import numpy as np\r
import matplotlib.cm as cm\r
import matplotlib.patches as patch\r
import matplotlib.colors as col\r
\r
-FOUT = 'dimer_in_potential.png'\r
-FCOF = 'rk_out.dat'\r
-FP1 = 'rk_point1.dat'\r
-FP2 = 'rk_point2.dat'\r
+FOUT = '../dimer_in_potential.png'\r
+FCOF = '../rk_out.dat'\r
+FP1 = '../rk_point1.dat'\r
+FP2 = '../rk_point2.dat'\r
\r
L = 2*np.pi\r
w = 0.1\r
tau = L/w\r
-evry = 10 # plot periodicity, don't set < 10 (bug)\r
-off = -10 # offset of (x,y) in L\r
+evry = 1 # plot periodicity, don't set < 10 (bug)\r
+offx = -98 # offset of (x,y) in L\r
+offy = 94\r
+tfrom = 90 # from t\r
+tto = 100 # to t\r
+driveangle = 41 # angle of ac drive\r
sper = 4 # spacial periodicity\r
pad = 0.01 # small padding so the ticks get drawn correctly...\r
msize = 20\r
show_plot = False\r
xlabel = r'$x$'\r
ylabel = r'$y$'\r
+cblabel = r'$\cos(x)\cos(y)+\cos(x)+\cos(y)$'\r
fs_ticks = 14\r
fs_labels = 24\r
title = 'F = 0.1, alpha = 83, l = 3.0, a = 1.5, eta1 = 1.0, eta2 = 1.5, theta = 0'\r
\r
print "reading data..."\r
\r
-data_cof = mlab.csv2rec(FCOF, delimiter='\t', comments='#')\r
+data_cof= mlab.csv2rec(FCOF, delimiter='\t', comments='#')\r
data_p1 = mlab.csv2rec(FP1, delimiter='\t', comments='#')\r
data_p2 = mlab.csv2rec(FP2, delimiter='\t', comments='#')\r
\r
\r
print "done."\r
\r
-tstep = t[1] - t[0]\r
+tstep = t[10] - t[9]\r
\r
-lim_min = off*L\r
-lim_max = (off+sper)*L\r
+fromt = tfrom/tstep\r
+tot = tto/tstep\r
\r
-xticks = np.arange(lim_min,lim_max+pad,L)\r
-yticks = np.arange(lim_min,lim_max+pad,L)\r
-xlim = (lim_min-pad,lim_max+pad)\r
-ylim = (lim_min-pad,lim_max+pad)\r
+every = evry*tstep\r
+trunc_cofx = cofx[fromt*tau:tot*tau:evry]\r
+trunc_cofy = cofy[fromt*tau:tot*tau:evry]\r
+trunc_p1x = p1x[fromt*tau:tot*tau:evry]\r
+trunc_p2x = p2x[fromt*tau:tot*tau:evry]\r
+trunc_p2y = p2y[fromt*tau:tot*tau:evry]\r
+trunc_p1y = p1y[fromt*tau:tot*tau:evry]\r
+trunc_pxdiff = pxdiff[fromt*tau:tot*tau:evry]\r
+trunc_pydiff = pydiff[fromt*tau:tot*tau:evry]\r
+trunc_t = t[fromt*tau:tot*tau:evry]\r
+\r
+xmax = np.ceil(max(trunc_cofx)/L)\r
+xmin = np.ceil(min(trunc_cofx)/L)\r
+ymax = np.ceil(max(trunc_cofy)/L)\r
+ymin = np.ceil(min(trunc_cofy)/L)\r
+\r
+print xmin*L, xmax*L\r
+xdiff = np.ceil(abs(max(trunc_cofx)-min(trunc_cofx))/L)\r
+ydiff = np.ceil(abs(max(trunc_cofy)-min(trunc_cofy))/L)\r
+\r
+xcorr = xdiff-sper\r
+ycorr = ydiff-sper\r
+\r
+print xdiff, ydiff\r
+print xdiff-(xdiff-sper)\r
+\r
+if xmax > 0:\r
+ limx_max = xmax*L\r
+ limx_min = (xmax-sper)*L\r
+if xmax < 0:\r
+ limx_min = xmin*L\r
+ limx_max = (xmin+sper)*L\r
+if ymax > 0:\r
+ limy_max = ymax*L\r
+ limy_min = (ymax-sper)*L\r
+if ymax < 0:\r
+ limy_min = ymin*L\r
+ limy_max = (ymin+sper)*L\r
+\r
+print (limx_max-limx_min)/L\r
+\r
+print "x:","(",limx_min,",",limx_max,")"\r
+print "y:","(",limy_min,",",limy_max,")"\r
+\r
+xticks = np.arange(limx_min,limx_max+pad,L)\r
+yticks = np.arange(limy_min,limy_max+pad,L)\r
+xlim = (limx_min-pad,limx_max+pad)\r
+ylim = (limy_min-pad,limy_max+pad)\r
tlabel = ("-2L","-L","0","L","2L")\r
\r
print "generating basemap..."\r
delta = 0.01\r
-x = np.arange(lim_min,lim_max, delta)\r
-y = np.arange(lim_min,lim_max, delta)\r
+x = np.arange(limx_min,limx_max, delta)\r
+y = np.arange(limy_min,limy_max, delta)\r
X,Y = np.meshgrid(x, y)\r
-Z = np.cos(X)*np.cos(Y)\r
+Z = np.cos(X)*np.cos(Y)+np.cos(X)+np.cos(Y)\r
print "done."\r
\r
fig = plt.figure(figsize=res)\r
rgb = ls.shade(Z,cmap)\r
print "done."\r
print "plotting basemap..."\r
- pot = plt.imshow(rgb, extent=[lim_min,lim_max,lim_min,lim_max])\r
+ pot = plt.imshow(rgb, extent=[limx_min,limx_max,limy_min,limy_max])\r
print "done."\r
else:\r
print "plotting basemap..."\r
print "done."\r
\r
cb = plt.colorbar(pot)\r
-cb.set_label(r'$\cos(x)\cos(y)$')\r
+cb.set_label(cblabel)\r
\r
-every = evry*tstep\r
-\r
-k = 0\r
print "plotting dimer onto the basemap..."\r
-cof = ax.scatter(cofx[::evry], cofy[::evry], c='yellow', marker='o', lw=0, s=msize/4, label='cof')\r
-x1 = ax.scatter(p1x[::evry], p1y[::evry], c='blue', marker='o', lw=0, s=msize, label='x1')\r
-x2 = ax.scatter(p2x[::evry], p2y[::evry], c='purple', marker='o', lw=0, s=msize, label='x2') \r
-for i in range(len(t)):\r
- if (t[i] == k*every and i <= len(t)):\r
- #print "plotting because",t[i],"=",k,"*",every,"(",k*every,")"\r
- bond = ax.arrow(p1x[i], p1y[i], pxdiff[i], pydiff[i], color='r', lw=0.1, label='bond')\r
- k = k+1\r
- else:\r
- #print "not plotting because",t[i],"!=",k,"*",every,"(",k*every,")"\r
- continue\r
+cof = ax.scatter(trunc_cofx, trunc_cofy, c='yellow', marker='o', lw=0, s=msize/4, label='cof')\r
+x1 = ax.scatter(trunc_p1x, trunc_p1y, c='blue', marker='o', lw=0, s=msize, label='x1')\r
+x2 = ax.scatter(trunc_p2x, trunc_p2y, c='purple', marker='o', lw=0, s=msize, label='x2') \r
+#bond = ax.arrow(trunc_p1x, trunc_p1y, trunc_pxdiff, trunc_pydiff, color='r', lw=0.1, label='bond')\r
+\r
+\r
+k=0\r
+for i in range(len(trunc_t)):\r
+ #bond = ax.arrow(p1x[i], p1y[i], pxdiff[i], pydiff[i], color='r', lw=0.1, label='bond')\r
+ if trunc_t[i] >= k*every:\r
+ #print "plotting because",t[i],"=",k,"*",every,"(",k*every,")"\r
+ bond = ax.arrow(trunc_p1x[i], trunc_p1y[i],\r
+ trunc_pxdiff[i], trunc_pydiff[i],\r
+ color='r', lw=0.1, label='bond',\r
+ alpha=0.5)\r
+ k = k+1\r
+ elif trunc_t[i] < k*every:\r
+ #print "not plotting because",t[i],"!=",k,"*",every,"(",k*every,")"\r
+ continue\r
+\r
print "done."\r
\r
+ax.arrow(limx_min,limy_min, 10*np.cos(driveangle*np.pi/180),\r
+ 10*np.sin(driveangle*np.pi/180), lw=5.0, color='g')\r
+\r
+"""\r
+n = 0\r
+for i in range(len(x)):\r
+ if x[i] >= n*L:\r
+ s1 = ax.axhline(y=n/2)\r
+ n = n+1\r
+ elif x[i] < n*L:\r
+ continue\r
+"""\r
+\r
+\r
ax.set_xlabel(xlabel, size=fs_labels)\r
ax.set_ylabel(ylabel, size=fs_labels)\r
ax.set_xticks(xticks)\r
ax.set_yticks(yticks)\r
-ax.set_xlim(xlim)\r
-ax.set_ylim(ylim)\r
+ax.set_xlim(limx_min,limx_max)\r
+ax.set_ylim(limy_min,limy_max)\r
ax.set_xticklabels(tlabel, size=fs_ticks)\r
ax.set_yticklabels(tlabel, size=fs_ticks)\r
\r
-legend = fig.legend((cof,x1,x2,bond),('cof','x1','x2','bond'), shadow=True, fancybox=True)\r
+legend = fig.legend((cof, x1, x2, bond), (r'$\vec{r}_s$', r'$\vec{r}_1$',\r
+ r'$\vec{r}_2$', r'$|\vec{r}_1-\vec{r}_2|$'), shadow=True,\r
+ fancybox=True)\r
legend.get_frame().set_alpha(0.75)\r
\r
-plt.title(title)\r
+#plt.title(title)\r
\r
print "saving to file.."\r
plt.savefig(FOUT)\r
print "showing plot..."\r
plt.show()\r
print "done."\r
+\r
diff --git a/quplot.h b/quplot.h
--- a/quplot.h
+++ b/quplot.h
* 38 = tor2
* 39 = stepr2
* 40 = theta
+ * 41 = psi
+ * 42 = U - amplitude of potential
*/
QString str1;