1 //=============================================================================
4 // contains all the magic math
6 //=============================================================================
9 //== INCLUDES =================================================================
11 #include "calculator.h"
14 //== IMPLEMENTATION ===========================================================
17 Calculator::Calculator(QuPlot
*parent
) :
25 Calculator::~Calculator()
30 //-----------------------------------------------------------------------------
31 void Calculator::derive_init()
33 file_rk_bif
= "../rk_bifurcate.dat";
34 file_rk_avg
= "../rk_bif_avg_velo.dat";
35 file_rk_out
= "../rk_out.dat";
36 file_rk_pt1
= "../rk_point1.dat";
37 file_rk_pt2
= "../rk_point2.dat";
38 file_rk_biv
= "../rk_bif_velo.dat";
40 file_eu_bif
= "../eu_bifurcate.dat";
41 file_eu_avg
= "../eu_bif_avg_velo.dat";
42 file_eu_out
= "../eu_out.dat";
43 file_eu_pt1
= "../eu_point1.dat";
44 file_eu_pt2
= "../eu_point2.dat";
45 file_eu_biv
= "../eu_bif_velo.dat";
47 file_rk_conv
= "../rk_conv.dat";
48 file_eu_conv
= "../eu_conv.dat";
50 file_rk_variance
= "../rk_variance.dat";
51 file_rk_distance
= "../rk_msd.dat";
52 v
= c
->valuesMap
.value(1);
56 void Calculator::derivs(const Doub x
, VecDoub_I
&y
, VecDoub_O
&dydx
,
57 const Doub params
[], Normaldev
&gauss
)
76 //static Doub G, xi, xi1, xi2, c, s, st, ax, ay, e0, e1, e2, e3, eta, l1, l2;
77 static Doub xi
, G
, s
, c
, st
, ax
, ay
;
78 static Doub x1
, y1
, x2
, y2
, F1x
, F1y
, F2x
, F2y
, fx
, fy
, stx
, sty
;
79 //for (Int i = 0;i<10;i++)
80 // qDebug() << "params["<<i<<"] :" << params[i];
85 G
= sqrt(2.0*params
[5]*params
[9])*xi
;
87 dydx
[1] = -params
[5]*y
[1]-sin(y
[0])+params
[1]*sin(params
[4]*x
)+params
[2]+G
;
91 dydx
[1] = -params
[5]*y
[1]-y
[0];
103 eta = params[5]+params[6];
105 e1 = 1.0/(params[5]*params[3]);
106 e2 = 1.0/(params[6]*params[3]);
107 e3 = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
109 l1 = params[3]*e0*params[5];
110 l2 = params[3]*e0*params[6];
114 //qDebug() << eta << "=" << params[5] << "+" << params[6];
115 //qDebug() << e0 << "=" << "1.0/"<< eta;
116 //qDebug() << e1 << "=" << "1.0/"<<params[5]<<"*"<<params[3];
117 //qDebug() << e2 << "=" << "1.0/"<<params[6]<<"*"<<params[3];
118 //qDebug() << e3 << "=" << "1.0/"<<params[5]<<"-"<<"1.0/"<<params[6];
120 //G = sqrt(2.0*eta*params[9]);
126 x1
= y
[0] + params
[7]*c
;
127 y1
= y
[1] + params
[7]*s
;
128 x2
= y
[0] - params
[8]*c
;
129 y2
= y
[1] - params
[8]*s
;
131 st
= sin(params
[4]*x
);
133 //qDebug() << st << "=" << "sin("<<params[4]<<"*"<<x<<")";
135 ax
= cos(params
[1])*params
[0]*st
;
136 ay
= sin(params
[1])*params
[0]*st
;
138 //fx = cos(params[10])*params[2];
139 //fy = sin(params[10])*params[2];
144 //qDebug() << "ax =" << ax << "=" << params[0] << "*" << "cos("<<params[1]<<")";
145 //qDebug() << "ay =" << ay << "=" << params[0] << "*" << "sin("<<params[1]<<")";
146 //qDebug() << "fx =" << fx << "=" << "cos("<<params[10]<<")*"<<params[2];
147 //qDebug() << "fy =" << fy << "=" << "sin("<<params[10]<<")*"<<params[2];
148 //qDebug() << "stx =" << stx <<"="<< ax << "*" << st << "+"<< fx;
149 //qDebug() << "sty =" << sty <<"="<< ay << "*" << st << "+" << fy;
153 //qDebug() << y[0] + params[7]*c << "=" << y[0] << "+" << params[7] << "*" << c;
154 //qDebug() << y[1] + params[7]*c << "=" << y[1] << "+" << params[7] << "*" << s;
155 //qDebug() << y[0] - params[8]*c << "=" << y[0] << "-" << params[8] << "*" << c;
156 //qDebug() << y[1] - params[8]*c << "=" << y[1] << "-" << params[8] << "*" << s;
158 F1x
= sin(x1
)*cos(y1
);//+1.6*sin(x1)+1.0*cos(y1);
159 F1y
= cos(x1
)*sin(y1
);//+1.6*cos(x1)+1.0*sin(y1);
160 F2x
= sin(x2
)*cos(y2
);//+1.6*sin(x2)+1.0*cos(y1);
161 F2y
= cos(x2
)*sin(y2
);//+1.6*cos(x2)+1.0*sin(y2);
163 dydx
[0] = (F1x
+F2x
+2.0*stx
)*params
[11];
164 dydx
[1] = (F1y
+F2y
+2.0*sty
)*params
[11];
165 dydx
[2] = params
[12]*(c
*F1y
-s
*F1x
)-params
[13]*(c
*F2y
-s
*F2x
)+params
[14]*(c
*sty
-s
*stx
);
171 eta = params[5]+params[6];
172 l1 = params[3]*params[5]/eta;
173 l2 = params[3]*params[6]/eta;
176 F1x = sin(x1)*cos(y1);//+1.6*sin(x1)+1.0*cos(y1);
177 F1y = cos(x1)*sin(y1);//+1.6*cos(x1)+1.0*sin(y1);
178 F2x = sin(x2)*cos(y2);//+1.6*sin(x2)+1.0*cos(y1);
179 F2y = cos(x2)*sin(y2);//+1.6*cos(x2)+1.0*sin(y2);
181 dydx[0] = (F1x+F2x+2.0*(params[0]*cos(params[1])*sin(params[4]*x)+params[10]*params[2])/eta;
182 dydx[1] = (F1y+F2y+2.0*(params[0]*sin(params[1])*sin(params[4]*x)+params[10]*params[2])/eta;
183 dydx[2] = (-s*F1x+c*F1y)/(params[5]*params[3]) - (-s*F2x+c*F2y)/(params[6]*params[3])
184 + (1.0/params[5] - 1.0/params[6])*params[0]*( (-s*cos(params[1])+c*sin(params[1])+sin(params[4]*x))+(- );
191 //-----------------------------------------------------------------------------
194 void Calculator::rk4(VecDoub_I
&y
, VecDoub_I
&dydx
, const Doub x
, const Doub h
,
195 VecDoub_O
&yout
, const Doub params
[], Normaldev
&gauss
)
198 VecDoub
dym(n
),dyt(n
),yt(n
);
202 for (Int i
=0;i
<n
;i
++) yt
[i
] = y
[i
]+hh
*dydx
[i
];
203 derivs(xh
,yt
,dyt
,params
,gauss
);
204 for (Int i
=0;i
<n
;i
++) yt
[i
] = y
[i
]+hh
*dyt
[i
];
205 derivs(xh
,yt
,dym
,params
,gauss
);
206 for (Int i
=0;i
<n
;i
++) {
207 yt
[i
] = y
[i
]+h
*dym
[i
];
210 derivs(x
+h
,yt
,dyt
,params
,gauss
);
211 for (Int i
=0;i
<n
;i
++) {
212 yout
[i
] = y
[i
]+h6
*(dydx
[i
]+dyt
[i
]+2.0*dym
[i
]);
217 void Calculator::euler(VecDoub_I
&y
, VecDoub_I
&dydx
, const Doub x
,
218 const Doub h
, VecDoub_O
&yout
)
221 for (Int i
=0;i
<n
;i
++) {
222 yout
[i
] = y
[i
]+h
*dydx
[i
];
225 //-----------------------------------------------------------------------------
229 * returns a random number between 0 and 2*pi
231 double Calculator::randpi(Ran
&r
)
233 Doub s
= 2.0*M_PI
*r
.doub();
237 //-----------------------------------------------------------------------------
240 * Writes the of two maps into a tab-separated file
241 * The second map *must* be specified.
243 void Calculator::write_to_file(const QString
& name
,
244 const QMap
<double, double> & map_x
,
245 const QMap
<double, double> & map_y
,
246 const QMap
<double, double> & map_phi
,
249 qDebug() << "writing to file" << name
;
253 QTextStream
out(&file
);
254 file
.open(QIODevice::WriteOnly
);
256 for (Int i
=0;i
<10;i
++)
257 out
<< "params["<<i
<<"] = " << params
[i
] << " ";
259 out
<< "r\tx\ty\tphi" << endl
;
260 QMapIterator
<double,double> x(map_x
);
261 QMapIterator
<double,double> y(map_y
);
262 QMapIterator
<double,double> phi(map_phi
);
265 while (x
.hasNext()) {
267 out
<< x
.key() << "\t" << x
.value() << endl
;
271 while (x
.hasNext() && y
.hasNext() && phi
.hasNext()) {
275 out
<< x
.key() << "\t" << x
.value() << "\t" << y
.value() << "\t" << phi
.value() << endl
;
280 qDebug () << "done.";
284 //-----------------------------------------------------------------------------
287 Doub
Calculator::mean(VecDoub_I
&msd
, const Int
& samples
)
290 VecDoub
var(samples
,0.0);
292 assert(msd
.size() == samples
);
294 for(Int i
=0;i
<samples
;i
++) {
299 for(Int i
=0;i
<samples
;i
++) {
300 var
[i
] = SQR(msd
[i
])-SQR(mu
);
302 assert(var
.size() == samples
) ;
304 for(Int j
=0;j
<samples
;j
++) {
315 * This function is kinda redundant right now, but will be needed later
317 void Calculator::gen_gnuplot(const QString
& name
)
320 qDebug() << "generating gnuplot file";
323 QTextStream
out(&file
);
324 file
.open(QIODevice::WriteOnly
);
328 out
<< "set lmargin 10;";
329 out
<< "set rmargin 3;";
330 out
<< "set multiplot;";
331 out
<< "set tmargin = 3;";
332 out
<< "set bmargin = 0;";
333 out
<< "set border 2+4+8;";
334 out
<< "set size 1,0.66667;";
335 out
<< "set origin 0.0,0.33333;";
337 out
<< "set format x \"\";";
338 out
<< "unset xlabel;";
339 out
<< "set ylabel \"x(kT) mod L \";";
340 out
<< "set ytics (\"0\" 0, \"L/2\", \"L\" L);";
341 out
<< "set yrange [0:L];";
342 out
<< "plot \""<<file_rk_bif
<<"\" w dots 1;";
343 out
<< "unset title;";
344 out
<< "set tmargin 0;";
345 out
<< "set bmargin 3;";
346 out
<< "unset format;";
347 out
<< "set size 1,0.33333;";
348 out
<< "set border 1+2+4+8;";
349 out
<< "set origin 0.0,0.0;";
350 out
<< "set xtics 0,0.1;";
351 out
<< "set xlabel \"a\";";
352 out
<< "set ylabel \"v\";";
353 out
<< "set yrange [-0.75:0.75];";
354 out
<< "set ytics -0.5,0.5;";
355 out
<< "plot \""<<file_rk_avg
<<"\" w dots 3;";
356 out
<< "unset multiplot;";
360 void Calculator::filter(QMap
<double, double> & map
, const Doub
& eps
)
362 qDebug() << "filter crap smaller than " << eps
;
365 QList
<double> values
, list
;
366 list
= map
.uniqueKeys();
369 //qDebug() << "size of map:" << size;
370 for(Int i
= 0; i
< size
; i
++) {
372 //qDebug() << "["<<key<<"]" << "writing values into list";
373 values
= map
.values(key
);
374 //qDebug() << "["<<key<<"]" << values.size() << "values for" << key;
376 for(Int k
=0;k
<values
.size()-1;k
++) {
377 diff
= abs(values
.at(k
+1) - values
.at(k
));
378 //qDebug() << values.at(k+1) <<"-"<< values.at(k) << "=" << diff;
380 //qDebug() << diff << "<=" << eps;
381 //qDebug() << "["<<k<<"]" << "removing" << values.at(k) << "from list";
382 values
.replace(k
+1,values
.at(k
));
385 for (int j
= 0; j
< values
.size(); j
++)
386 //qDebug() << "["<<j<<"]" << values.at(j);
387 //qDebug() << "removing" << key << "from map";
389 for (int j
= 0; j
< values
.size(); j
++) {
390 //qDebug() << "inserting" << values.at(j) << "into map at" << key;
391 map
.insertMulti(key
,values
.at(j
));
393 //qDebug() << map.values(key);
394 //qDebug() << endl << "end while" << endl;
398 void Calculator::IC_setup(VecDoub
&y
, Doub x0
, Doub xdot0
, Doub y0
, Doub ydot0
,
399 Doub phi0
, Doub params
[])
402 const Int sys
= c
->valuesMap
.value(1);
405 qDebug() << "using " << "\tx0 = " << x0 << "\ty0 = " << y0
406 << "\txdot0 = " << xdot0 << "\tydot0 = " << ydot0
407 << "\tphi0 = " << fmod(phi0,360) << endl;
412 params
[0] = c
->valuesMap
.value(2); // a;
413 params
[1] = c
->valuesMap
.value(3); // A;
414 params
[2] = c
->valuesMap
.value(4); // F;
415 params
[3] = c
->valuesMap
.value(5); // l;
416 params
[4] = c
->valuesMap
.value(6); // w;
417 params
[5] = c
->valuesMap
.value(25); // eta1;
419 params
[9] = c
->valuesMap
.value(32); // T;
421 y
[0] = x0
, y
[1] = xdot0
;
424 params
[0] = c
->valuesMap
.value(2); // a;
425 params
[1] = fmod(c
->valuesMap
.value(3),360)*PI
/180; // alpha;
426 params
[2] = c
->valuesMap
.value(4); // F;
427 params
[3] = c
->valuesMap
.value(5); // l;
428 params
[4] = c
->valuesMap
.value(6); // w;
429 params
[5] = c
->valuesMap
.value(25); // eta1;
430 params
[6] = c
->valuesMap
.value(26); // eta2;
431 params
[9] = c
->valuesMap
.value(32); // T;
432 params
[10] = fmod(c
->valuesMap
.value(40),360)*PI
/180; // theta;
433 Doub e0
= 1.0/(params
[5]+params
[6]);
434 Doub e1
= 1.0/(params
[5]*params
[3]);
435 Doub e2
= 1.0/(params
[6]*params
[3]);
436 Doub e3
= 1.0/params
[3]*((1.0/params
[5])-(1.0/params
[6]));
437 Doub fx
= cos(params
[10])*params
[2];
438 Doub fy
= sin(params
[10])*params
[2];
439 params
[11] = e0
; // 1/eta
445 params
[7] = params
[3]*params
[11]*params
[5]; // l1;
446 params
[8] = params
[3]*params
[11]*params
[6]; // l2;
447 y
[0] = x0
, y
[1] = y0
, y
[2] = phi0
;
451 for (Int i = 0; i<n;i++)
452 qDebug() << "y("<<i<<") =" << y[i];
453 qDebug() << "eta1 = " << eta1 << "\teta2 = " << eta2 << "\tA = " << A << "\tF = " << F;
457 void Calculator::bifparam(const int &v
, const QString
¶m
, const Doub
&r
,
464 else if (param
== "F")
466 else if (param
== "A")
469 throw ("Fatal error!");
473 params
[1] = fmod(r
,360)*PI
/180;
474 else if (param
== "F") {
476 params
[15] = cos(params
[10])*params
[2];
477 params
[16] = sin(params
[10])*params
[2];
479 else if (param
== "a")
481 else if (param
== "w")
483 else if (param
== "l") {
485 params
[11] = 1.0/(params
[5]+params
[6]);
486 params
[12] = 1.0/(params
[5]*params
[3]);
487 params
[13] = 1.0/(params
[6]*params
[3]);
488 params
[14] = 1.0/params
[3]*((1.0/params
[5])-(1.0/params
[6]));
489 params
[7] = params
[3]*params
[11]*params
[5]; // l1;
490 params
[8] = params
[3]*params
[11]*params
[6]; // l2;
492 else if (param
== "e") {
494 params
[11] = 1.0/(params
[5]+params
[6]);
495 params
[12] = 1.0/(params
[5]*params
[3]);
496 params
[13] = 1.0/(params
[6]*params
[3]);
497 params
[14] = 1.0/params
[3]*((1.0/params
[5])-(1.0/params
[6]));
498 params
[7] = params
[3]*params
[11]*params
[5]; // l1;
499 params
[8] = params
[3]*params
[11]*params
[6]; // l2;
501 else if (param
== "t")
502 params
[10] = fmod(r
,360)*PI
/180;
504 throw("Fatal error!");
506 qDebug() << "Bifurcation parameter is " << param
;
510 void Calculator::set_color(const Doub vel
, QString
& color
)
513 if (vel
> 1.0-ctol
&& vel
< 1.0+ctol
)
515 else if (vel
> 0.5-ctol
&& vel
< 0.5+ctol
)
517 else if (vel
> 0.0-ctol
&& vel
< 0.0+ctol
)
519 else if (vel
> -0.5-ctol
&& vel
< -0.5+ctol
)
521 else if (vel
> -1.0-ctol
&& vel
< -1.0+ctol
)
523 else color
= "black";
526 void Calculator::calc2(const Int n
)
528 // setting up a few start values
529 qDebug() << "calc2 has been called: " << "calc2("<<n
<<")";
531 Normaldev
gauss(0,1,time(NULL
));
533 Doub x1
,x2
,y1
,y2
,xtest
;
535 Doub x0
, xdot0
, y0
, ydot0
, phi0
;
536 x0
= c
->valuesMap
.value(8);
537 xdot0
= c
->valuesMap
.value(10);
538 ydot0
= c
->valuesMap
.value(24);
539 y0
= c
->valuesMap
.value(23);
540 phi0
= c
->valuesMap
.value(27);
541 const Int samples
= c
->valuesMap
.value(35);
542 const Int v
= c
->valuesMap
.value(1);
543 const Int M
= c
->valuesMap
.value(11);
544 const Doub fromt
= c
->valuesMap
.value(22);
545 const Int s
= c
->valuesMap
.value(14);
546 const Doub h
= c
->valuesMap
.value(13);
547 const Doub t0
= c
->valuesMap
.value(7);
548 const Doub fromr1
= c
->valuesMap
.value(17);
549 const Doub fromr2
= c
->valuesMap
.value(37);
550 const Doub tor1
= c
->valuesMap
.value(18);
551 const Doub tor2
= c
->valuesMap
.value(38);
552 const Doub stepr1
= c
->valuesMap
.value(19);
553 const Doub stepr2
= c
->valuesMap
.value(39);
554 const Doub eps
= c
->valuesMap
.value(36);
555 const Doub from
= c
->valuesMap
.value(15);
556 const Doub to
= c
->valuesMap
.value(16);
559 const bool ruk
= c
->ruk
;
560 const bool eul
= c
->eul
;
561 QString param1
= c
->str1
;
562 QString param2
= c
->str2
;
563 qDebug() << "System" << v
;
567 Doub xpi
[samples
][s
];
568 Doub ypi
[samples
][s
];
569 Doub phipi
[samples
][s
];
578 qDebug() << "standard mode";
579 qDebug() << "stepsize:" << h
;
580 qDebug() << "number of steps:" << s
;
582 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
585 tau
= 2*PI
/params
[4];
588 if (xtest
>= per
*tau
) {
591 derivs(x
,y
,dydx
,params
,gauss
); // sets first dydx[] for rk4()
592 rk4(y
,dydx
,x
,h
,yout
,params
,gauss
);
593 if (per
>= s
*fromt
) {
596 rkValues_x
.insert(x
,y
[0]);
599 //xpi = ABS(fmod(yout[0],(2*PI)));
600 //ypi = ABS(fmod(yout[1],(2*PI)));
601 rkValues_x
.insert(x
,y
[0]);
602 rkValues_y
.insert(x
,y
[1]);
603 //rkValues_x.insert(x,xpi);
604 //rkValues_y.insert(x,ypi);
605 rkValues_phi
.insert(x
,y
[2]);
606 x1
= y
[0] + params
[7]*cos(y
[2]);
607 y1
= y
[1] + params
[7]*sin(y
[2]);
608 x2
= y
[0] - params
[8]*cos(y
[2]);
609 y2
= y
[1] - params
[8]*sin(y
[2]);
610 rkValues_p1x
.insert(x
,x1
);
611 rkValues_p1y
.insert(x
,y1
);
612 rkValues_p2x
.insert(x
,x2
);
613 rkValues_p2y
.insert(x
,y2
);
617 for (Int j
=0;j
<n
;j
++)
621 qDebug() << "sizeof(rkValues_x) =" << rkValues_x
.size();
622 qDebug() << "sizeof(rkValues_y) =" << rkValues_y
.size();
623 qDebug() << "sizeof(rkValues_phi) =" << rkValues_phi
.size();
624 write_to_file(file_rk_out
, rkValues_x
, rkValues_y
, rkValues_phi
, params
);
625 write_to_file(file_rk_pt1
, rkValues_p1x
, rkValues_p1y
, rkValues_phi
, params
);
626 write_to_file(file_rk_pt2
, rkValues_p2x
, rkValues_p2y
, rkValues_phi
, params
);
629 qDebug() << "converge mode";
632 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
633 for (iter
=1;(iter
*from
)<to
;iter
++) {
637 qDebug() << "stepsize is" << hh
<< "with" << ss
<< "total steps.";
638 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
639 for (Int i
=0;i
<ss
&&x
<(ss
-0.5);i
++) {
640 derivs(x
,y
,dydx
,params
,gauss
);
641 rk4(y
,dydx
,x
,hh
,yout
,params
,gauss
);
642 for (Int j
=0;j
<n
;j
++)
648 rkValues_x
.insert(hh
,yout
[0]);
651 rkValues_x
.insert(hh
,yout
[0]);
652 rkValues_y
.insert(hh
,yout
[1]);
653 rkValues_phi
.insert(hh
,yout
[2]);
657 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
658 for (Int i
=0;i
<ss
&&x
<(ss
-0.5);i
++) {
659 derivs(x
,y
,dydx
,params
,gauss
);
660 euler(y
,dydx
,x
,hh
,yout
);
661 for (Int j
=0;j
<n
;j
++)
667 euValues_x
.insert(hh
,yout
[0]);
670 euValues_x
.insert(hh
,yout
[0]);
671 euValues_y
.insert(hh
,yout
[1]);
672 euValues_phi
.insert(hh
,yout
[2]);
676 qDebug() << "sizeof(rkValues_x) =" << rkValues_x
.size();
677 write_to_file(file_rk_conv
, rkValues_x
, rkValues_y
, rkValues_phi
, params
);
678 write_to_file(file_eu_conv
, euValues_x
, euValues_y
, euValues_phi
, params
);
681 qDebug() << "bifurcation mode";
682 qDebug() << "bifurcation parameter: " << param1
;
683 qDebug() << "Number of samples : " << samples
;
685 for (Int j
=0;(j
*stepr1
)<=(tor1
-fromr1
);j
++) {
687 qDebug() << "at step" << j
<< "of" << (tor1
-fromr1
)/stepr1
;
688 for (Int z
=0;z
<samples
;z
++) {
690 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
691 //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
692 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
693 bifparam(v
,param1
,r
,params
);
694 //qDebug() << x << y[0] << y[1] << y[2];
695 //for(Int i = 0;i<=10;i++)
696 // qDebug() << params[i];
697 // inner for: iterate over runge-kutta
701 tau
= 2*PI
/params
[4];
702 //qDebug() << "tau: " << tau << "\tL: " << L;
704 xtest
= x
+h
; // test next x increment
705 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
706 if (xtest
>= k
*tau
) {
707 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
708 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
710 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
711 derivs(x
,y
,dydx
,params
,gauss
);
712 rk4(y
,dydx
,x
,hnew
,yout
,params
,gauss
);
715 xpi
[z
][k
] = ABS(fmod(yout
[0],(2*PI
)));
716 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
718 rkValues_x
.insert(r
,xpi
[z
][k
]);
719 rkValues_xdot
.insert(r
,dydx
[0]);
723 xpi
[z
][k
] = ABS(fmod(yout
[0],(2*PI
)));
724 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
725 ypi
[z
][k
] = ABS(fmod(yout
[1],(2*PI
)));
726 phipi
[z
][k
] = ABS(fmod(yout
[2],(2*PI
)));
727 rkDist_x
.insert(r
,SQR(yout
[0]-x0
));
728 rkDist_y
.insert(r
,SQR(yout
[1]-y0
));
729 rkDist_phi
.insert(r
,SQR(yout
[2]-phi0
));
731 rkValues_x
.insert(r
,xpi
[z
][k
]);
732 rkValues_y
.insert(r
,ypi
[z
][k
]);
733 rkValues_phi
.insert(r
,yout
[2]);
734 rkValues_xdot
.insert(r
,dydx
[0]);
735 rkValues_ydot
.insert(r
,dydx
[1]);
736 rkValues_phidot
.insert(r
,dydx
[2]);
742 derivs(x
,y
,dydx
,params
,gauss
);
743 rk4(y
,dydx
,x
,h
,yout
,params
,gauss
); // xold
744 for (Int j
=0;j
<n
;j
++)
745 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
746 x
+= h
; // do it for real
748 //qDebug() << "at" << x << "in time for k = " << k;
753 avg
[0] = (Doub
) m_x
/k
;
754 rkAvgvelo_x
.insert(r
,avg
[0]);
757 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
758 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
; avg
[2] = (Doub
) m_phi
/k
;
759 rkAvgvelo_x
.insert(r
,avg
[0]);
760 rkAvgvelo_y
.insert(r
,avg
[1]);
761 rkAvgvelo_phi
.insert(r
,avg
[2]);
765 VecDoub
msd_x(samples
,0.0);
766 VecDoub
msd_y(samples
,0.0);
767 VecDoub
msd_phi(samples
,0.0);
768 VecDoub
sigma_x(s
,0.0);
769 VecDoub
sigma_y(s
,0.0);
770 VecDoub
sigma_phi(s
,0.0);
771 //qDebug() << "samples:" << samples << "steps:" << s;
772 for(Int i
=0;i
<s
;i
++) {
773 //qDebug() << "["<<i<<"]:-----------------------------------------";
774 for(Int j
=0;j
<samples
;j
++) {
775 //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
776 msd_x
[j
] = xpi
[j
][i
];
777 //qDebug() << "["<<j<<"]:" << sper_x[j];
778 msd_y
[j
] = ypi
[j
][i
];
779 msd_phi
[j
] = phipi
[j
][i
];
781 sigma_x
[i
] = mean(msd_x
,samples
);
782 sigma_y
[i
] = mean(msd_y
,samples
);
783 sigma_phi
[i
] = mean(msd_phi
,samples
);
784 //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
786 for(Int k
=0;k
<s
;k
++) {
787 rkVariance_x
.insert(r
,sigma_x
[k
]);
788 rkVariance_y
.insert(r
,sigma_y
[k
]);
789 rkVariance_phi
.insert(r
,sigma_phi
[k
]);
791 //qDebug() << rkVariance_x;
794 * we need to do something about this
797 filter(rkValues_x,eps);
798 filter(rkValues_y,eps);
799 filter(rkValues_xdot,eps);
800 filter(rkValues_ydot,eps);
801 filter(rkAvgvelo_x, eps);
802 filter(rkAvgvelo_y, eps);
803 filter(rkAvgvelo_phi,eps);*/
804 qDebug() << "sizeof(rkValues_x) =" << rkValues_x
.size();
805 qDebug() << "sizeof(rkValues_y) =" << rkValues_y
.size();
806 qDebug() << "sizeof(rkValues_xdot) =" << rkValues_xdot
.size();
807 qDebug() << "sizeof(rkValues_ydot) =" << rkValues_ydot
.size();
808 write_to_file(file_rk_bif
, rkValues_x
, rkValues_y
, rkValues_phi
, params
);
809 write_to_file(file_rk_biv
, rkValues_xdot
, rkValues_ydot
, rkValues_phidot
, params
);
810 write_to_file(file_rk_avg
, rkAvgvelo_x
, rkAvgvelo_y
, rkAvgvelo_phi
, params
);
811 write_to_file(file_rk_variance
, rkVariance_x
, rkVariance_y
, rkVariance_phi
,params
);
812 write_to_file(file_rk_distance
, rkDist_x
, rkDist_y
, rkDist_phi
, params
);
815 qDebug() << "dual bifurcation mode";
816 qDebug() << "bifurcation parameters: " << param1
<<"\t"<<param2
;
817 qDebug() << "Number of samples : " << samples
;
819 QString vx_color
, vy_color
;
820 QString name
= "dual_bif.dat";
822 QTextStream
out(&file
);
823 file
.open(QIODevice::WriteOnly
);
824 out
<< "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl
;
825 for(Int it1
=0;(it1
*stepr1
)<=(tor1
-fromr1
);it1
++) {
826 for(Int it2
=0;(it2
*stepr2
)<=(tor2
-fromr2
);it2
++) {
827 stp1
= fromr1
+it1
*stepr1
;
828 stp2
= fromr2
+it2
*stepr2
;
829 qDebug() << stp1
<< "("<<tor1
<<")" << stp2
<< "("<<tor2
<<")";
830 for (Int z
=0;z
<samples
;z
++) {
832 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
833 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
834 bifparam(v
,param1
,stp1
,params
);
835 bifparam(v
,param2
,stp2
,params
);
838 for (Int i=0;i<10;i++)
839 out << "params["<<i<<"] = " << params[i] << " ";
842 // inner for: iterate over runge-kutta
846 tau
= 2*PI
/params
[4];
847 //qDebug() << "tau: " << tau << "\tL: " << L;
849 xtest
= x
+h
; // test next x increment
850 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
851 if (xtest
>= k
*tau
) {
852 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
853 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
855 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
856 derivs(x
,y
,dydx
,params
,gauss
);
857 rk4(y
,dydx
,x
,hnew
,yout
,params
,gauss
);
860 derivs(x
,y
,dydx
,params
,gauss
);
861 rk4(y
,dydx
,x
,h
,yout
,params
,gauss
); // xold
862 for (Int j
=0;j
<n
;j
++)
863 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
864 x
+= h
; // do it for real
867 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
868 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
;
869 avg
[2] = (Doub
) m_phi
/k
;
870 set_color(avg
[0], vx_color
);
871 set_color(avg
[1], vy_color
);
872 out
<< stp1
<< "\t" << stp2
<< "\t"
873 << avg
[0] << "\t" << vx_color
<< "\t" << avg
[1]
874 << "\t" << vy_color
<< "\t" << avg
[2] << endl
;
876 } // end 2. outer for
877 } // end 1. outer for
886 qDebug() << "(euler) standard mode";
887 qDebug() << "stepsize:" << h
;
888 qDebug() << "number of steps:" << s
;
890 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
893 tau
= 2*PI
/params
[4];
896 if (xtest
>= per
*tau
) {
899 derivs(x
,y
,dydx
,params
,gauss
); // sets first dydx[] for rk4()
900 euler(y
,dydx
,x
,h
,yout
);
901 if (per
>= s
*fromt
) {
904 euValues_x
.insert(x
,y
[0]);
907 //xpi = ABS(fmod(yout[0],(2*PI)));
908 //ypi = ABS(fmod(yout[1],(2*PI)));
909 euValues_x
.insert(x
,y
[0]);
910 euValues_y
.insert(x
,y
[1]);
911 //rkValues_x.insert(x,xpi);
912 //rkValues_y.insert(x,ypi);
913 euValues_phi
.insert(x
,y
[2]);
914 x1
= y
[0] + params
[7]*cos(y
[2]);
915 y1
= y
[1] + params
[7]*sin(y
[2]);
916 x2
= y
[0] - params
[8]*cos(y
[2]);
917 y2
= y
[1] - params
[8]*sin(y
[2]);
918 euValues_p1x
.insert(x
,x1
);
919 euValues_p1y
.insert(x
,y1
);
920 euValues_p2x
.insert(x
,x2
);
921 euValues_p2y
.insert(x
,y2
);
925 for (Int j
=0;j
<n
;j
++)
929 qDebug() << "sizeof(euValues_x) =" << euValues_x
.size();
930 qDebug() << "sizeof(euValues_y) =" << euValues_y
.size();
931 qDebug() << "sizeof(euValues_phi) =" << euValues_phi
.size();
932 write_to_file(file_eu_out
, euValues_x
, euValues_y
, euValues_phi
, params
);
933 write_to_file(file_eu_pt1
, euValues_p1x
, euValues_p1y
, euValues_phi
, params
);
934 write_to_file(file_eu_pt2
, euValues_p2x
, euValues_p2y
, euValues_phi
, params
);
937 qDebug() << "bifurcation mode";
938 qDebug() << "bifurcation parameter: " << param1
;
939 qDebug() << "Number of samples : " << samples
;
941 for (Int j
=0;(j
*stepr1
)<=(tor1
-fromr1
);j
++) {
943 qDebug() << "at step" << j
<< "of" << (tor1
-fromr1
)/stepr1
;
944 for (Int z
=0;z
<samples
;z
++) {
946 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
947 //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
948 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
949 bifparam(v
,param1
,r
,params
);
950 // inner for: iterate over runge-kutta
954 tau
= 2*PI
/params
[4];
955 //qDebug() << "tau: " << tau << "\tL: " << L;
957 xtest
= x
+h
; // test next x increment
958 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
959 if (xtest
>= k
*tau
) {
960 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
961 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
963 //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
964 derivs(x
,y
,dydx
,params
,gauss
);
965 euler(y
,dydx
,x
,hnew
,yout
);
969 xpi
[z
][k
] = ABS(fmod(yout
[0],(2*PI
)));
970 euValues_x
.insert(r
,xpi
[z
][k
]);
971 euValues_xdot
.insert(r
,dydx
[0]);
974 xpi
[z
][k
] = ABS(fmod(yout
[0],(2*PI
)));
975 ypi
[z
][k
] = ABS(fmod(yout
[1],(2*PI
)));
976 euValues_x
.insert(r
,xpi
[z
][k
]);
977 euValues_y
.insert(r
,ypi
[z
][k
]);
978 euValues_phi
.insert(r
,yout
[2]);
979 euValues_xdot
.insert(r
,dydx
[0]);
980 euValues_ydot
.insert(r
,dydx
[1]);
981 euValues_phidot
.insert(r
,dydx
[2]);
987 derivs(x
,y
,dydx
,params
,gauss
);
988 euler(y
,dydx
,x
,h
,yout
); // xold
989 for (Int j
=0;j
<n
;j
++)
990 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
991 x
+= h
; // do it for real
993 //qDebug() << "at" << x << "in time for k = " << k;
998 avg
[0] = (Doub
) m_x
/k
;
999 euAvgvelo_x
.insert(r
,avg
[0]);
1002 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
1003 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
; avg
[2] = (Doub
) m_phi
/k
;
1004 euAvgvelo_x
.insert(r
,avg
[0]);
1005 euAvgvelo_y
.insert(r
,avg
[1]);
1006 euAvgvelo_phi
.insert(r
,avg
[2]);
1012 * we need to do something about this
1014 filter(euValues_x
,eps
);
1015 filter(euValues_y
,eps
);
1016 filter(euValues_xdot
,eps
);
1017 filter(euValues_ydot
,eps
);
1018 filter(euAvgvelo_x
, eps
);
1019 filter(euAvgvelo_y
, eps
);
1020 filter(euAvgvelo_phi
,eps
);
1021 qDebug() << "sizeof(euValues_x) =" << euValues_x
.size();
1022 qDebug() << "sizeof(euValues_y) =" << euValues_y
.size();
1023 qDebug() << "sizeof(euValues_xdot) =" << euValues_xdot
.size();
1024 qDebug() << "sizeof(euValues_ydot) =" << euValues_ydot
.size();
1025 write_to_file(file_eu_bif
, euValues_x
, euValues_y
, euValues_phi
, params
);
1026 write_to_file(file_eu_biv
, euValues_xdot
, euValues_ydot
, euValues_phidot
, params
);
1027 write_to_file(file_eu_avg
, euAvgvelo_x
, euAvgvelo_y
, euAvgvelo_phi
, params
);
1030 qDebug() << "dual bifurcation mode";
1031 qDebug() << "bifurcation parameters: " << param1
<<"\t"<<param2
;
1032 qDebug() << "Number of samples : " << samples
;
1034 QString name
= "dual_bif.dat";
1036 QTextStream
out(&file
);
1037 file
.open(QIODevice::WriteOnly
);
1038 for (Int it1
=0;(it1
*stepr1
)<=(tor1
-fromr1
);it1
++) {
1039 for (Int it2
=0;(it2
*stepr2
)<=(tor2
-fromr1
);it2
++) {
1040 stp1
= fromr1
+it1
*stepr1
;
1041 stp2
= fromr2
+it2
*stepr2
;
1042 qDebug() << stp1
<< "("<<tor1
<<")" << stp2
<< "("<<tor2
<<")";
1043 for (Int z
=0;z
<samples
;z
++) {
1045 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
1046 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
1047 bifparam(v
,param1
,stp1
,params
);
1048 bifparam(v
,param2
,stp2
,params
);
1049 // inner for: iterate over runge-kutta
1053 tau
= 2*PI
/params
[4];
1054 //qDebug() << "tau: " << tau << "\tL: " << L;
1056 xtest
= x
+h
; // test next x increment
1057 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1058 if (xtest
>= k
*tau
) {
1059 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1060 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1062 //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1063 derivs(x
,y
,dydx
,params
,gauss
);
1064 euler(y
,dydx
,x
,hnew
,yout
);
1067 derivs(x
,y
,dydx
,params
,gauss
);
1068 euler(y
,dydx
,x
,h
,yout
); // xold
1069 for (Int j
=0;j
<n
;j
++)
1070 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
1071 x
+= h
; // do it for real
1073 Int m_x
, m_y
, m_phi
;
1074 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
1075 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
; avg
[2] = (Doub
) m_phi
/k
;
1076 out
<< scientific
<< stp1
<< "\t" << stp2
<< "\t" << avg
[0] << "\t" << avg
[1] << "\t" << avg
[2] << endl
;
1078 } // end 2. outer for
1079 } // end 1. outer for
1086 throw("No method specified!");
1089 rkValues_xdot
.clear();
1091 rkValues_ydot
.clear();
1092 rkValues_phi
.clear();
1093 rkValues_phidot
.clear();
1094 rkAvgvelo_x
.clear();
1095 rkAvgvelo_y
.clear();
1096 rkAvgvelo_phi
.clear();
1099 euValues_xdot
.clear();
1101 euValues_ydot
.clear();
1102 euValues_phi
.clear();
1103 euValues_phidot
.clear();
1104 euAvgvelo_x
.clear();
1105 euAvgvelo_y
.clear();
1106 euAvgvelo_phi
.clear();
1110 //=============================================================================