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 file_rk_defangle
= "../rk_defangle.dat";
53 file_eu_variance
= "../eu_variance.dat";
54 file_eu_distance
= "../eu_msd.dat";
55 file_eu_defangle
= "../eu_defangle.dat";
57 file_rk_dualbif
= "../rk_dualbif.dat";
58 file_eu_dualbif
= "../eu_dualbif.dat";
60 v
= c
->valuesMap
.value(1);
64 void Calculator::derivs(const Doub x
, VecDoub_I
&y
, VecDoub_O
&dydx
,
86 //static Doub G, xi, xi1, xi2, c, s, st, ax, ay, e0, e1, e2, e3, eta, l1, l2;
87 static Doub xi
, G
, s
, c
, st
, ax
, ay
, cp
, sp
, m1
, m2
, p1
, p2
;
88 static Doub x1
, y1
, x2
, y2
, F1x
, F1y
, F2x
, F2y
, stx
, sty
;
89 //for (Int i = 0;i<10;i++)
90 // qDebug() << "params["<<i<<"] :" << params[i];
94 G
= sqrt(2.0*params
[5]*params
[9])*xi
;
96 dydx
[1] = -params
[5]*y
[1]-sin(y
[0])+params
[1]*sin(params
[4]*x
)+params
[2]+G
;
100 dydx
[1] = -params
[5]*y
[1]-y
[0];
115 eta = params[5]+params[6];
117 e1 = 1.0/(params[5]*params[3]);
118 e2 = 1.0/(params[6]*params[3]);
119 e3 = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
121 l1 = params[3]*e0*params[5];
122 l2 = params[3]*e0*params[6];
126 //qDebug() << eta << "=" << params[5] << "+" << params[6];
127 //qDebug() << e0 << "=" << "1.0/"<< eta;
128 //qDebug() << e1 << "=" << "1.0/"<<params[5]<<"*"<<params[3];
129 //qDebug() << e2 << "=" << "1.0/"<<params[6]<<"*"<<params[3];
130 //qDebug() << e3 << "=" << "1.0/"<<params[5]<<"-"<<"1.0/"<<params[6];
132 //G = sqrt(2.0*eta*params[9]);
138 x1
= y
[0] + params
[7]*c
;
139 y1
= y
[1] + params
[7]*s
;
140 x2
= y
[0] - params
[8]*c
;
141 y2
= y
[1] - params
[8]*s
;
143 st
= sin(params
[4]*x
);
145 //qDebug() << st << "=" << "sin("<<params[4]<<"*"<<x<<")";
147 ax
= cos(params
[1])*params
[0]*st
;
148 ay
= sin(params
[1])*params
[0]*st
;
150 //fx = cos(params[10])*params[2];
151 //fy = sin(params[10])*params[2];
156 //qDebug() << "ax =" << ax << "=" << params[0] << "*" << "cos("<<params[1]<<")";
157 //qDebug() << "ay =" << ay << "=" << params[0] << "*" << "sin("<<params[1]<<")";
158 //qDebug() << "fx =" << fx << "=" << "cos("<<params[10]<<")*"<<params[2];
159 //qDebug() << "fy =" << fy << "=" << "sin("<<params[10]<<")*"<<params[2];
160 //qDebug() << "stx =" << stx <<"="<< ax << "*" << st << "+"<< fx;
161 //qDebug() << "sty =" << sty <<"="<< ay << "*" << st << "+" << fy;
165 //qDebug() << y[0] + params[7]*c << "=" << y[0] << "+" << params[7] << "*" << c;
166 //qDebug() << y[1] + params[7]*c << "=" << y[1] << "+" << params[7] << "*" << s;
167 //qDebug() << y[0] - params[8]*c << "=" << y[0] << "-" << params[8] << "*" << c;
168 //qDebug() << y[1] - params[8]*c << "=" << y[1] << "-" << params[8] << "*" << s;
170 cp
= cos(params
[17]);
171 sp
= sin(params
[17]);
178 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);
179 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);
180 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);
181 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);
183 dydx
[0] = (F1x
+F2x
+2.0*stx
/*+sqrt(2*params[9])*(xix1+xix2)*/)*params
[11];
184 dydx
[1] = (F1y
+F2y
+2.0*sty
/*+sqrt(2*params[9])*(xiy1+xiy2)*/)*params
[11];
185 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))*/;
191 void Calculator::derivw(const Doub x
, VecDoub_I
&y
, VecDoub_O
&dW
,
192 const Doub params
[], Normaldev
&gauss
)
204 //params[10] = theta;
214 Doub xix1
, xix2
, xiy1
, xiy2
;
215 Doub c
, s
, T
, l
, eta1
, eta2
;
223 l
= params
[7] + params
[8];
232 dW
[0] = sqrt(2*T
)*(sqrt(eta1
)*xix1
+sqrt(eta2
)*xix2
);
233 dW
[1] = sqrt(2*T
)*(sqrt(eta1
)*xiy1
+sqrt(eta2
)*xiy2
);
234 dW
[2] = (sqrt(2*T
)/l
)*((1/sqrt(eta1
))*(-s
*xix1
+c
*xiy1
)-(1/sqrt(eta2
))*(-s
*xix2
+c
*xiy2
));
238 //-----------------------------------------------------------------------------
241 void Calculator::rk4(VecDoub_I
&y
, VecDoub_I
&dydx
, const Doub x
, const Doub h
,
242 VecDoub_O
&yout
, const Doub params
[])
245 VecDoub
dym(n
),dyt(n
),yt(n
);
249 for (Int i
=0;i
<n
;i
++) yt
[i
] = y
[i
]+hh
*dydx
[i
];
250 derivs(xh
,yt
,dyt
,params
);
251 for (Int i
=0;i
<n
;i
++) yt
[i
] = y
[i
]+hh
*dyt
[i
];
252 derivs(xh
,yt
,dym
,params
);
253 for (Int i
=0;i
<n
;i
++) {
254 yt
[i
] = y
[i
]+h
*dym
[i
];
257 derivs(x
+h
,yt
,dyt
,params
);
258 for (Int i
=0;i
<n
;i
++) {
259 yout
[i
] = y
[i
]+h6
*(dydx
[i
]+dyt
[i
]+2.0*dym
[i
]);
264 void Calculator::euler(VecDoub_I
&y
, VecDoub_I
&dydx
, VecDoub_I
&dW
, const Doub x
,
265 const Doub h
, VecDoub_O
&yout
)
268 for (Int i
=0;i
<n
;i
++) {
269 yout
[i
] = y
[i
]+h
*dydx
[i
]+sqrt(h
)*dW
[i
];
272 //-----------------------------------------------------------------------------
276 * returns a random number between 0 and 2*pi
278 double Calculator::randpi(Ran
&r
)
280 Doub s
= 2.0*M_PI
*r
.doub();
284 //-----------------------------------------------------------------------------
287 * Writes the of two maps into a tab-separated file
288 * The second map *must* be specified.
290 void Calculator::write_to_file(const QString
& name
,
291 const QMap
<double, double> & map_x
,
292 const QMap
<double, double> & map_y
,
293 const QMap
<double, double> & map_phi
,
296 qDebug() << "writing to file" << name
;
300 QTextStream
out(&file
);
301 file
.open(QIODevice::WriteOnly
);
303 for (Int i
=0;i
<10;i
++)
304 out
<< "params["<<i
<<"] = " << params
[i
] << " ";
306 out
<< "r\tx\ty\tphi" << endl
;
307 QMapIterator
<double,double> x(map_x
);
308 QMapIterator
<double,double> y(map_y
);
309 QMapIterator
<double,double> phi(map_phi
);
312 while (x
.hasNext()) {
314 out
<< x
.key() << "\t" << x
.value() << endl
;
318 if (map_x
.size() > 0 && map_y
.size() == 0 && map_phi
.size() == 0) {
319 while (x
.hasNext()) {
321 out
<< x
.key() << "\t" << x
.value() << endl
;
324 else if (map_x
.size() > 0 && map_y
.size() > 0 && map_phi
.size() > 0) {
325 while (x
.hasNext() && y
.hasNext() && phi
.hasNext()) {
329 //qDebug() << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value();
330 out
<< x
.key() << "\t" << x
.value() << "\t" << y
.value() << "\t" << phi
.value() << endl
;
336 qDebug () << "done.";
340 //-----------------------------------------------------------------------------
343 void Calculator::diffusion(VecDoub_I
&msd
, const Int
& samples
, Doub
& sigma
)
346 VecDoub
var(samples
,0.0);
349 for(Int i
=0;i
<samples
;i
++) {
350 //qDebug() << "msd["<<i<<"]" << msd[i];
354 //qDebug() << "mu:" << mu;
356 for(Int i
=0;i
<samples
;i
++) {
357 var
[i
] = SQR(msd
[i
])-SQR(mu
);
358 //qDebug() << "var["<<i<<"]" << var[i];
361 for(Int j
=0;j
<samples
;j
++) {
367 //qDebug() << "sigma:" << sigma;
371 * This function is kinda redundant right now, but will be needed later
373 void Calculator::gen_gnuplot(const QString
& name
)
376 qDebug() << "generating gnuplot file";
379 QTextStream
out(&file
);
380 file
.open(QIODevice::WriteOnly
);
384 out
<< "set lmargin 10;";
385 out
<< "set rmargin 3;";
386 out
<< "set multiplot;";
387 out
<< "set tmargin = 3;";
388 out
<< "set bmargin = 0;";
389 out
<< "set border 2+4+8;";
390 out
<< "set size 1,0.66667;";
391 out
<< "set origin 0.0,0.33333;";
393 out
<< "set format x \"\";";
394 out
<< "unset xlabel;";
395 out
<< "set ylabel \"x(kT) mod L \";";
396 out
<< "set ytics (\"0\" 0, \"L/2\", \"L\" L);";
397 out
<< "set yrange [0:L];";
398 out
<< "plot \""<<file_rk_bif
<<"\" w dots 1;";
399 out
<< "unset title;";
400 out
<< "set tmargin 0;";
401 out
<< "set bmargin 3;";
402 out
<< "unset format;";
403 out
<< "set size 1,0.33333;";
404 out
<< "set border 1+2+4+8;";
405 out
<< "set origin 0.0,0.0;";
406 out
<< "set xtics 0,0.1;";
407 out
<< "set xlabel \"a\";";
408 out
<< "set ylabel \"v\";";
409 out
<< "set yrange [-0.75:0.75];";
410 out
<< "set ytics -0.5,0.5;";
411 out
<< "plot \""<<file_rk_avg
<<"\" w dots 3;";
412 out
<< "unset multiplot;";
416 void Calculator::filter(QMap
<double, double> & map
, const Doub
& eps
)
418 qDebug() << "filter crap smaller than " << eps
;
421 QList
<double> values
, list
;
422 list
= map
.uniqueKeys();
425 //qDebug() << "size of map:" << size;
426 for(Int i
= 0; i
< size
; i
++) {
428 //qDebug() << "["<<key<<"]" << "writing values into list";
429 values
= map
.values(key
);
430 //qDebug() << "["<<key<<"]" << values.size() << "values for" << key;
432 for(Int k
=0;k
<values
.size()-1;k
++) {
433 diff
= abs(values
.at(k
+1) - values
.at(k
));
434 //qDebug() << values.at(k+1) <<"-"<< values.at(k) << "=" << diff;
436 //qDebug() << diff << "<=" << eps;
437 //qDebug() << "["<<k<<"]" << "removing" << values.at(k) << "from list";
438 values
.replace(k
+1,values
.at(k
));
441 for (int j
= 0; j
< values
.size(); j
++)
442 //qDebug() << "["<<j<<"]" << values.at(j);
443 //qDebug() << "removing" << key << "from map";
445 for (int j
= 0; j
< values
.size(); j
++) {
446 //qDebug() << "inserting" << values.at(j) << "into map at" << key;
447 map
.insertMulti(key
,values
.at(j
));
449 //qDebug() << map.values(key);
450 //qDebug() << endl << "end while" << endl;
454 void Calculator::IC_setup(VecDoub
&y
, Doub x0
, Doub y0
, Doub xdot0
, Doub ydot0
,
455 Doub phi0
, Doub params
[])
458 const Int sys
= c
->valuesMap
.value(1);
461 qDebug() << "using " << "\tx0 = " << x0 << "\ty0 = " << y0
462 << "\txdot0 = " << xdot0 << "\tydot0 = " << ydot0
463 << "\tphi0 = " << fmod(phi0,360) << endl;
468 params
[0] = c
->valuesMap
.value(2); // a;
469 params
[1] = c
->valuesMap
.value(3); // A;
470 params
[2] = c
->valuesMap
.value(4); // F;
471 params
[3] = c
->valuesMap
.value(5); // l;
472 params
[4] = c
->valuesMap
.value(6); // w;
473 params
[5] = c
->valuesMap
.value(25); // eta1;
475 params
[9] = c
->valuesMap
.value(32); // T;
477 y
[0] = x0
, y
[1] = xdot0
;
480 params
[0] = c
->valuesMap
.value(2); // a;
481 params
[1] = fmod(c
->valuesMap
.value(3),360)*PI
/180; // alpha;
482 params
[2] = c
->valuesMap
.value(4); // F;
483 params
[3] = c
->valuesMap
.value(5); // l;
484 params
[4] = c
->valuesMap
.value(6); // w;
485 params
[5] = c
->valuesMap
.value(25); // eta1;
486 params
[6] = c
->valuesMap
.value(26); // eta2;
487 params
[9] = c
->valuesMap
.value(32); // T;
488 params
[10] = fmod(c
->valuesMap
.value(40),360)*PI
/180; // theta;
489 params
[17] = fmod(c
->valuesMap
.value(41),360)*PI
/180; // psi;
490 params
[18] = c
->valuesMap
.value(42); // U;
491 Doub e0
= 1.0/(params
[5]+params
[6]);
492 Doub e1
= 1.0/(params
[5]*params
[3]);
493 Doub e2
= 1.0/(params
[6]*params
[3]);
494 Doub e3
= 1.0/params
[3]*((1.0/params
[5])-(1.0/params
[6]));
495 Doub fx
= cos(params
[10])*params
[2];
496 Doub fy
= sin(params
[10])*params
[2];
497 params
[11] = e0
; // 1/eta
503 params
[7] = params
[3]*params
[11]*params
[5]; // l1;
504 params
[8] = params
[3]*params
[11]*params
[6]; // l2;
505 y
[0] = x0
, y
[1] = y0
, y
[2] = phi0
;
509 for (Int i = 0; i<n;i++)
510 qDebug() << "y("<<i<<") =" << y[i];
511 qDebug() << "eta1 = " << eta1 << "\teta2 = " << eta2 << "\tA = " << A << "\tF = " << F;
515 void Calculator::bifparam(const int &v
, const QString
¶m
, const Doub
&r
,
522 else if (param
== "F")
524 else if (param
== "A")
527 throw ("Fatal error!");
531 params
[1] = fmod(r
,360)*PI
/180;
532 else if (param
== "F") {
534 params
[15] = cos(params
[10])*params
[2];
535 params
[16] = sin(params
[10])*params
[2];
537 else if (param
== "a")
539 else if (param
== "w")
541 else if (param
== "l") {
543 params
[11] = 1.0/(params
[5]+params
[6]);
544 params
[12] = 1.0/(params
[5]*params
[3]);
545 params
[13] = 1.0/(params
[6]*params
[3]);
546 params
[14] = 1.0/params
[3]*((1.0/params
[5])-(1.0/params
[6]));
547 params
[7] = params
[3]*params
[11]*params
[5]; // l1;
548 params
[8] = params
[3]*params
[11]*params
[6]; // l2;
550 else if (param
== "e") {
552 params
[11] = 1.0/(params
[5]+params
[6]);
553 params
[12] = 1.0/(params
[5]*params
[3]);
554 params
[13] = 1.0/(params
[6]*params
[3]);
555 params
[14] = 1.0/params
[3]*((1.0/params
[5])-(1.0/params
[6]));
556 params
[7] = params
[3]*params
[11]*params
[5]; // l1;
557 params
[8] = params
[3]*params
[11]*params
[6]; // l2;
559 else if (param
== "t") {
560 params
[10] = fmod(r
,360)*PI
/180;
561 params
[15] = cos(params
[10])*params
[2];
562 params
[16] = sin(params
[10])*params
[2];
565 throw("Fatal error!");
567 qDebug() << "Bifurcation parameter is " << param
;
571 void Calculator::set_color(const Doub vel
, QString
& color
)
574 if (vel
> 1.0-ctol
&& vel
< 1.0+ctol
)
576 else if (vel
> 0.5-ctol
&& vel
< 0.5+ctol
)
578 else if (vel
> 0.0-ctol
&& vel
< 0.0+ctol
)
580 else if (vel
> -0.5-ctol
&& vel
< -0.5+ctol
)
582 else if (vel
> -1.0-ctol
&& vel
< -1.0+ctol
)
584 else color
= "black";
587 void Calculator::calc2(const Int n
)
589 // setting up a few start values
590 qDebug() << "calc2 has been called: " << "calc2("<<n
<<")";
591 Normaldev
gauss(0,1,time(NULL
));
593 Doub x1
,x2
,y1
,y2
,xtest
;
595 Doub x0
, xdot0
, y0
, ydot0
, phi0
;
596 x0
= c
->valuesMap
.value(8);
597 xdot0
= c
->valuesMap
.value(10);
598 ydot0
= c
->valuesMap
.value(24);
599 y0
= c
->valuesMap
.value(23);
600 phi0
= c
->valuesMap
.value(27);
601 const Int samples
= c
->valuesMap
.value(35);
602 const Int v
= c
->valuesMap
.value(1);
603 const Int M
= c
->valuesMap
.value(11);
604 const Doub fromt
= c
->valuesMap
.value(22);
605 const Int s
= c
->valuesMap
.value(14);
606 const Doub h
= c
->valuesMap
.value(13);
607 const Doub t0
= c
->valuesMap
.value(7);
608 const Doub fromr1
= c
->valuesMap
.value(17);
609 const Doub fromr2
= c
->valuesMap
.value(37);
610 const Doub tor1
= c
->valuesMap
.value(18);
611 const Doub tor2
= c
->valuesMap
.value(38);
612 const Doub stepr1
= c
->valuesMap
.value(19);
613 const Doub stepr2
= c
->valuesMap
.value(39);
614 const Doub eps
= c
->valuesMap
.value(36);
615 const Doub from
= c
->valuesMap
.value(15);
616 const Doub to
= c
->valuesMap
.value(16);
617 const Doub hstep
= c
->valuesMap
.value(43);
620 const bool ruk
= c
->ruk
;
621 const bool eul
= c
->eul
;
622 QString param1
= c
->str1
;
623 QString param2
= c
->str2
;
624 qDebug() << "System" << v
;
627 Doub trajectory_length
, deflection
, velocity_angle
, bias_angle
;
630 Doub xpi
[samples
][s
];
631 Doub ypi
[samples
][s
];
632 Doub phipi
[samples
][s
];
642 qDebug() << "standard mode";
643 qDebug() << "stepsize:" << h
;
644 qDebug() << "number of steps:" << s
;
646 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
649 tau
= 2*PI
/params
[4];
652 if (xtest
>= k
*tau
) {
655 derivs(x
,y
,dydx
,params
); // sets first dydx[] for rk4()
656 rk4(y
,dydx
,x
,h
,yout
,params
);
657 for (Int j
=0;j
<n
;j
++)
663 rkValues_x
.insert(x
,yout
[0]);
666 //xpi = ABS(fmod(yout[0],(2*PI)));
667 //ypi = ABS(fmod(yout[1],(2*PI)));
668 rkValues_x
.insert(x
,yout
[0]);
669 rkValues_y
.insert(x
,yout
[1]);
670 //rkValues_x.insert(x,xpi);
671 //rkValues_y.insert(x,ypi);
672 rkValues_phi
.insert(x
,(180/PI
)*ABS(fmod(yout
[2],(2*PI
))));
673 x1
= y
[0] + params
[7]*cos(yout
[2]);
674 y1
= y
[1] + params
[7]*sin(yout
[2]);
675 x2
= y
[0] - params
[8]*cos(yout
[2]);
676 y2
= y
[1] - params
[8]*sin(yout
[2]);
677 rkValues_p1x
.insert(x
,x1
);
678 rkValues_p1y
.insert(x
,y1
);
679 rkValues_p2x
.insert(x
,x2
);
680 rkValues_p2y
.insert(x
,y2
);
685 qDebug() << "sizeof(rkValues_x) =" << rkValues_x
.size();
686 qDebug() << "sizeof(rkValues_y) =" << rkValues_y
.size();
687 qDebug() << "sizeof(rkValues_phi) =" << rkValues_phi
.size();
688 write_to_file(file_rk_out
, rkValues_x
, rkValues_y
, rkValues_phi
, params
);
689 write_to_file(file_rk_pt1
, rkValues_p1x
, rkValues_p1y
, rkValues_phi
, params
);
690 write_to_file(file_rk_pt2
, rkValues_p2x
, rkValues_p2y
, rkValues_phi
, params
);
693 qDebug() << "converge mode";
696 for (iter
=1;from
+(iter
*hstep
)<=to
;iter
++) {
697 for (Int z
=0;z
<samples
;z
++) {
698 qDebug() << "at step" << z
<< "of" << samples
;
700 hh
= from
+(iter
*hstep
);
702 qDebug() << "stepsize is" << hh
<< "with" << ss
<< "total steps.";
703 //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
704 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
707 tau
= 2*PI
/params
[4];
709 //qDebug() << k << "<" << s;
710 //qDebug() << k << tau << k*tau;
712 if (xtest
>= k
*tau
) {
713 //qDebug() << xtest << ">=" << k << "*" << tau << "("<<k*tau<<")";
715 derivs(x
,y
,dydx
,params
);
716 rk4(y
,dydx
,x
,hnew
,yout
,params
);
721 //qDebug() << xtest << "=" << x<<"+"<<hh;
722 derivs(x
,y
,dydx
,params
);
723 rk4(y
,dydx
,x
,hh
,yout
,params
); // xold
726 for (Int j
=0;j
<n
;j
++)
731 rkValues_x
.insert(hh
,yout
[0]);
734 rkValues_x
.insert(hh
,yout
[0]);
735 rkValues_y
.insert(hh
,yout
[1]);
736 rkValues_phi
.insert(hh
,yout
[2]);
740 //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
741 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
744 tau
= 2*PI
/params
[4];
747 if (xtest
>= k
*tau
) {
749 derivs(x
,y
,dydx
,params
);
750 derivw(x
,y
,dW
,params
,gauss
);
751 euler(y
,dydx
,dW
,x
,hnew
,yout
);
755 derivs(x
,y
,dydx
,params
);
756 derivw(x
,y
,dW
,params
,gauss
);
757 euler(y
,dydx
,dW
,x
,hh
,yout
);
760 for (Int j
=0;j
<n
;j
++)
765 euValues_x
.insert(hh
,yout
[0]);
768 euValues_x
.insert(hh
,yout
[0]);
769 euValues_y
.insert(hh
,yout
[1]);
770 euValues_phi
.insert(hh
,yout
[2]);
775 qDebug() << "sizeof(rkValues_x) =" << rkValues_x
.size();
776 write_to_file(file_rk_conv
, rkValues_x
, rkValues_y
, rkValues_phi
, params
);
777 write_to_file(file_eu_conv
, euValues_x
, euValues_y
, euValues_phi
, params
);
780 qDebug() << "bifurcation mode";
781 qDebug() << "bifurcation parameter: " << param1
;
782 qDebug() << "Number of samples : " << samples
;
784 for (Int j
=0;(j
*stepr1
)<=(tor1
-fromr1
);j
++) {
786 qDebug() << "at step" << j
<< "of" << (tor1
-fromr1
)/stepr1
;
787 for (Int z
=0;z
<samples
;z
++) {
789 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
790 //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
791 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
792 bifparam(v
,param1
,r
,params
);
793 //qDebug() << x << y[0] << y[1] << y[2];
794 //for(Int i = 0;i<=10;i++)
795 // qDebug() << params[i];
796 // inner for: iterate over runge-kutta
799 tau
= 2*PI
/params
[4];
801 //qDebug() << "tau: " << tau << "\tL: " << L;
803 xtest
= x
+h
; // test next x increment
804 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
805 if (xtest
>= k
*tau
) {
806 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
807 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
809 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
810 derivs(x
,y
,dydx
,params
);
811 rk4(y
,dydx
,x
,hnew
,yout
,params
);
814 xpi
[z
][k
] = ABS(fmod(yout
[0],(2*PI
)));
815 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
817 rkValues_x
.insert(r
,xpi
[z
][k
]);
818 rkValues_xdot
.insert(r
,dydx
[0]);
823 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
825 phipi
[z
][k
] = yout
[2];
826 rkDist_x
.insert(r
,SQR(yout
[0]-x0
));
827 rkDist_y
.insert(r
,SQR(yout
[1]-y0
));
828 rkDist_phi
.insert(r
,SQR(yout
[2]-phi0
));
830 rkValues_x
.insert(r
,ABS(fmod(xpi
[z
][k
],L
)));
831 rkValues_y
.insert(r
,ABS(fmod(ypi
[z
][k
],L
)));
832 rkValues_phi
.insert(r
,ABS(fmod(phipi
[z
][k
],L
)));
833 rkValues_xdot
.insert(r
,dydx
[0]);
834 rkValues_ydot
.insert(r
,dydx
[1]);
835 rkValues_phidot
.insert(r
,dydx
[2]);
841 derivs(x
,y
,dydx
,params
);
842 rk4(y
,dydx
,x
,h
,yout
,params
); // xold
843 for (Int j
=0;j
<n
;j
++)
844 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
845 x
+= h
; // do it for real
847 //qDebug() << "at" << x << "in time for k = " << k;
852 avg
[0] = (Doub
) m_x
/k
;
853 rkAvgvelo_x
.insert(r
,avg
[0]);
856 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
857 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
; avg
[2] = (Doub
) m_phi
/k
;
858 // makes only sense with r = params[10]
859 trajectory_length
= sqrt(SQR(yout
[0]-x0
)+SQR(yout
[1])-y0
);
860 if (abs(yout
[0]-x0
) > trajectory_length
) {
863 velocity_angle
= (acos((yout
[0]-x0
)/trajectory_length
));
865 bias_angle
= params
[10];
866 deflection
= (180/PI
)*ABS(velocity_angle
- bias_angle
);
868 qDebug() << velocity_angle << "("<<yout[0]-x0<<")"
869 << "-" << bias_angle << "="
871 << "[ length:"<<trajectory_length<<"]";
873 rkAvgvelo_x
.insert(r
,avg
[0]);
874 rkAvgvelo_y
.insert(r
,avg
[1]);
875 rkAvgvelo_phi
.insert(r
,avg
[2]);
876 rkDefAngle
.insert(r
,deflection
);
880 VecDoub
msd_x(samples
,0.0);
881 VecDoub
msd_y(samples
,0.0);
882 VecDoub
msd_phi(samples
,0.0);
886 //qDebug() << "samples:" << samples << "steps:" << s;
887 for(Int i
=0;i
<s
;i
++) {
888 //qDebug() << "["<<i<<"]:---------------------------------";
889 for(Int j
=0;j
<samples
;j
++) {
890 //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
891 msd_x
[j
] = xpi
[j
][i
];
892 //qDebug() << "["<<j<<"]:" << sper_x[j];
893 msd_y
[j
] = ypi
[j
][i
];
894 msd_phi
[j
] = phipi
[j
][i
];
896 diffusion(msd_x
, samples
, sigma
);
898 diffusion(msd_y
, samples
, sigma
);
900 diffusion(msd_phi
, samples
, sigma
);
901 sigma_phi
[i
] = sigma
;
902 //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
904 // the lists of keys in QMultiMap are read out backwards...
905 for(Int k
=(s
-1);k
>=0;k
--) {
906 rkVariance_x
.insert(r
,sigma_x
[k
]);
907 rkVariance_y
.insert(r
,sigma_y
[k
]);
908 rkVariance_phi
.insert(r
,sigma_phi
[k
]);
910 //qDebug() << rkVariance_x;
913 * we need to do something about this
916 filter(rkValues_x,eps);
917 filter(rkValues_y,eps);
918 filter(rkValues_xdot,eps);
919 filter(rkValues_ydot,eps);
920 filter(rkAvgvelo_x, eps);
921 filter(rkAvgvelo_y, eps);
922 filter(rkAvgvelo_phi,eps);*/
923 qDebug() << "sizeof(rkValues_x) =" << rkValues_x
.size();
924 qDebug() << "sizeof(rkValues_y) =" << rkValues_y
.size();
925 qDebug() << "sizeof(rkValues_xdot) =" << rkValues_xdot
.size();
926 qDebug() << "sizeof(rkValues_ydot) =" << rkValues_ydot
.size();
927 write_to_file(file_rk_bif
, rkValues_x
, rkValues_y
, rkValues_phi
, params
);
928 write_to_file(file_rk_biv
, rkValues_xdot
, rkValues_ydot
, rkValues_phidot
, params
);
929 write_to_file(file_rk_avg
, rkAvgvelo_x
, rkAvgvelo_y
, rkAvgvelo_phi
, params
);
930 write_to_file(file_rk_variance
, rkVariance_x
, rkVariance_y
, rkVariance_phi
,params
);
931 write_to_file(file_rk_distance
, rkDist_x
, rkDist_y
, rkDist_phi
, params
);
932 write_to_file(file_rk_defangle
, rkDefAngle
, NullMap
, NullMap
, params
);
935 qDebug() << "dual bifurcation mode";
936 qDebug() << "bifurcation parameters: " << param1
<<"\t"<<param2
;
937 qDebug() << "Number of samples : " << samples
;
939 QString vx_color
, vy_color
;
940 QString name
= "tmp.dat";
942 QTextStream
out(&file
);
943 file
.open(QIODevice::WriteOnly
);
944 out
<< "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl
;
945 for(Int it1
=0;(it1
*stepr1
)<=(tor1
-fromr1
);it1
++) {
946 for(Int it2
=0;(it2
*stepr2
)<=(tor2
-fromr2
);it2
++) {
947 stp1
= fromr1
+it1
*stepr1
;
948 stp2
= fromr2
+it2
*stepr2
;
949 qDebug() << stp1
<< "("<<tor1
<<")" << stp2
<< "("<<tor2
<<")";
950 for (Int z
=0;z
<samples
;z
++) {
952 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
953 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
954 bifparam(v
,param1
,stp1
,params
);
955 bifparam(v
,param2
,stp2
,params
);
958 for (Int i=0;i<10;i++)
959 out << "params["<<i<<"] = " << params[i] << " ";
962 // inner for: iterate over runge-kutta
965 tau
= 2*PI
/params
[4];
967 //qDebug() << "tau: " << tau << "\tL: " << L;
969 xtest
= x
+h
; // test next x increment
970 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
971 if (xtest
>= k
*tau
) {
972 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
973 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
975 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
976 derivs(x
,y
,dydx
,params
);
977 rk4(y
,dydx
,x
,hnew
,yout
,params
);
980 derivs(x
,y
,dydx
,params
);
981 rk4(y
,dydx
,x
,h
,yout
,params
); // xold
982 for (Int j
=0;j
<n
;j
++)
983 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
984 x
+= h
; // do it for real
987 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
988 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
;
989 avg
[2] = (Doub
) m_phi
/k
;
990 set_color(avg
[0], vx_color
);
991 set_color(avg
[1], vy_color
);
992 out
<< stp1
<< "\t" << stp2
<< "\t"
993 << avg
[0] << "\t" << vx_color
<< "\t" << avg
[1]
994 << "\t" << vy_color
<< "\t" << avg
[2] << endl
;
996 } // end 2. outer for
997 } // end 1. outer for
999 qDebug() << "moving" << name
<< "to safety";
1000 file
.copy(name
, file_rk_dualbif
);
1008 qDebug() << "(euler) standard mode";
1009 qDebug() << "stepsize:" << h
;
1010 qDebug() << "number of steps:" << s
;
1012 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
1016 tau
= 2*PI
/params
[4];
1018 //qDebug() << k << "<" << s;
1019 //qDebug() << "x is now" << x;
1020 //qDebug() << "deriving values at" << x;
1021 //qDebug() << x0, y0, phi0;
1022 derivs(x
,y
,dydx
,params
); // sets first dydx[] for rk4()
1023 derivw(x
,y
,dW
,params
,gauss
);
1024 euler(y
,dydx
,dW
,x
,h
,yout
);
1025 //for (Int i=0;i<=18;i++)
1026 //qDebug() << "params["<<i<<"]" << params[i];
1028 //qDebug() << x <<"=="<< hit*0.1;
1029 //qDebug() << "writing2 to map";
1031 //qDebug() << "next value to write out" << hit*0.1;
1035 euValues_x
.insert(x
,yout
[0]);
1038 //xpi = ABS(fmod(yout[0],(2*PI)));
1039 //ypi = ABS(fmod(yout[1],(2*PI)));
1040 euValues_x
.insert(x
,yout
[0]);
1041 euValues_y
.insert(x
,yout
[1]);
1042 //rkValues_x.insert(x,xpi);
1043 //rkValues_y.insert(x,ypi);
1044 euValues_phi
.insert(x
,yout
[2]);
1045 x1
= yout
[0] + params
[7]*cos(yout
[2]);
1046 y1
= yout
[1] + params
[7]*sin(yout
[2]);
1047 x2
= yout
[0] - params
[8]*cos(yout
[2]);
1048 y2
= yout
[1] - params
[8]*sin(yout
[2]);
1049 euValues_p1x
.insert(x
,x1
);
1050 euValues_p1y
.insert(x
,y1
);
1051 euValues_p2x
.insert(x
,x2
);
1052 euValues_p2y
.insert(x
,y2
);
1054 } // end fromt switch
1057 //qDebug() << "making step of size" << h;
1058 for (Int j
=0;j
<n
;j
++)
1061 //qDebug() << "new x is now" << x;
1063 //qDebug() << x << ">=" << k*tau;
1064 //qDebug() << "incrementing k";
1066 //qDebug() << "new value of x needed is" << k*tau;
1069 qDebug() << "sizeof(euValues_x) =" << euValues_x
.size();
1070 qDebug() << "sizeof(euValues_y) =" << euValues_y
.size();
1071 qDebug() << "sizeof(euValues_phi) =" << euValues_phi
.size();
1072 write_to_file(file_eu_out
, euValues_x
, euValues_y
, euValues_phi
, params
);
1073 write_to_file(file_eu_pt1
, euValues_p1x
, euValues_p1y
, euValues_phi
, params
);
1074 write_to_file(file_eu_pt2
, euValues_p2x
, euValues_p2y
, euValues_phi
, params
);
1077 qDebug() << "bifurcation mode";
1078 qDebug() << "bifurcation parameter: " << param1
;
1079 qDebug() << "Number of samples : " << samples
;
1081 for (Int j
=0;(j
*stepr1
)<=(tor1
-fromr1
);j
++) {
1082 r
= fromr1
+j
*stepr1
;
1083 qDebug() << "at step" << j
<< "of" << (tor1
-fromr1
)/stepr1
;
1084 for (Int z
=0;z
<samples
;z
++) {
1086 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
1087 //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
1088 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
1089 bifparam(v
,param1
,r
,params
);
1090 // inner for: iterate over runge-kutta
1093 tau
= 2*PI
/params
[4];
1095 //qDebug() << "tau: " << tau << "\tL: " << L;
1097 xtest
= x
+h
; // test next x increment
1098 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1099 if (xtest
>= k
*tau
) {
1100 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1101 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1103 //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1104 derivs(x
,y
,dydx
,params
);
1105 derivw(x
,y
, dW
, params
,gauss
);
1106 euler(y
,dydx
,dW
,x
,hnew
,yout
);
1109 xpi
[z
][k
] = yout
[0];
1111 euValues_x
.insert(r
,ABS(fmod(xpi
[z
][k
],L
)));
1112 euValues_xdot
.insert(r
,dydx
[0]);
1116 xpi
[z
][k
] = yout
[0];
1117 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
1118 ypi
[z
][k
] = yout
[1];
1119 phipi
[z
][k
] = yout
[2];
1120 euDist_x
.insert(r
,SQR(yout
[0]-x0
));
1121 euDist_y
.insert(r
,SQR(yout
[1]-y0
));
1122 euDist_phi
.insert(r
,SQR(yout
[2]-phi0
));
1124 euValues_x
.insert(r
,ABS(fmod(xpi
[z
][k
],L
)));
1125 euValues_y
.insert(r
,ABS(fmod(ypi
[z
][k
],L
)));
1126 euValues_phi
.insert(r
,ABS(fmod(phipi
[z
][k
],L
)));
1127 euValues_xdot
.insert(r
,dydx
[0]);
1128 euValues_ydot
.insert(r
,dydx
[1]);
1129 euValues_phidot
.insert(r
,dydx
[2]);
1135 derivs(x
,y
,dydx
,params
);
1136 derivw(x
,y
,dW
, params
,gauss
);
1137 euler(y
,dydx
,dW
,x
,h
,yout
); // xold
1138 for (Int j
=0;j
<n
;j
++)
1139 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
1140 x
+= h
; // do it for real
1142 //qDebug() << "at" << x << "in time for k = " << k;
1144 Int m_x
, m_y
, m_phi
;
1147 avg
[0] = (Doub
) m_x
/k
;
1148 euAvgvelo_x
.insert(r
,avg
[0]);
1151 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
1152 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
; avg
[2] = (Doub
) m_phi
/k
;
1153 trajectory_length
= sqrt(SQR(yout
[0]-x0
)+SQR(yout
[1])-y0
);
1154 if (abs(yout
[0]-x0
) > trajectory_length
) {
1155 velocity_angle
= PI
;
1157 velocity_angle
= (acos((yout
[0]-x0
)/trajectory_length
));
1159 bias_angle
= params
[10];
1160 deflection
= (180/PI
)*ABS(velocity_angle
- bias_angle
);
1161 euAvgvelo_x
.insert(r
,avg
[0]);
1162 euAvgvelo_y
.insert(r
,avg
[1]);
1163 euAvgvelo_phi
.insert(r
,avg
[2]);
1164 euDefAngle
.insert(r
,deflection
);
1168 VecDoub
msd_x(samples
,0.0);
1169 VecDoub
msd_y(samples
,0.0);
1170 VecDoub
msd_phi(samples
,0.0);
1173 Doub sigma_phi
[s
-1];
1174 //qDebug() << "samples:" << samples << "steps:" << s;
1175 for(Int i
=0;i
<s
;i
++) {
1176 //qDebug() << "["<<i<<"]:---------------------------------";
1177 for(Int j
=0;j
<samples
;j
++) {
1178 //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
1179 msd_x
[j
] = xpi
[j
][i
];
1180 //qDebug() << "["<<j<<"]:" << sper_x[j];
1181 msd_y
[j
] = ypi
[j
][i
];
1182 msd_phi
[j
] = phipi
[j
][i
];
1184 diffusion(msd_x
, samples
, sigma
);
1186 diffusion(msd_y
, samples
, sigma
);
1188 diffusion(msd_phi
, samples
, sigma
);
1189 sigma_phi
[i
] = sigma
;
1190 //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
1192 // the lists of keys in QMultiMap are read out backwards...
1193 for(Int k
=(s
-1);k
>=0;k
--) {
1194 euVariance_x
.insert(r
,sigma_x
[k
]);
1195 euVariance_y
.insert(r
,sigma_y
[k
]);
1196 euVariance_phi
.insert(r
,sigma_phi
[k
]);
1198 //qDebug() << euVariance_x;
1201 * we need to do something about this
1203 qDebug() << "sizeof(euValues_x) =" << euValues_x
.size();
1204 qDebug() << "sizeof(euValues_y) =" << euValues_y
.size();
1205 qDebug() << "sizeof(euValues_xdot) =" << euValues_xdot
.size();
1206 qDebug() << "sizeof(euValues_ydot) =" << euValues_ydot
.size();
1207 write_to_file(file_eu_bif
, euValues_x
, euValues_y
, euValues_phi
, params
);
1208 write_to_file(file_eu_biv
, euValues_xdot
, euValues_ydot
, euValues_phidot
, params
);
1209 write_to_file(file_eu_avg
, euAvgvelo_x
, euAvgvelo_y
, euAvgvelo_phi
, params
);
1210 write_to_file(file_eu_variance
, euVariance_x
, euVariance_y
, euVariance_phi
,params
);
1211 write_to_file(file_eu_distance
, euDist_x
, euDist_y
, euDist_phi
, params
);
1212 write_to_file(file_eu_defangle
, euDefAngle
, NullMap
, NullMap
, params
);
1215 qDebug() << "dual bifurcation mode";
1216 qDebug() << "bifurcation parameters: " << param1
<<"\t"<<param2
;
1217 qDebug() << "Number of samples : " << samples
;
1219 QString vx_color
, vy_color
;
1220 QString name
= "tmp.dat";
1222 QTextStream
out(&file
);
1223 file
.open(QIODevice::WriteOnly
);
1224 out
<< "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl
;
1225 for (Int it1
=0;(it1
*stepr1
)<=(tor1
-fromr1
);it1
++) {
1226 for (Int it2
=0;(it2
*stepr2
)<=(tor2
-fromr2
);it2
++) {
1227 stp1
= fromr1
+it1
*stepr1
;
1228 stp2
= fromr2
+it2
*stepr2
;
1229 qDebug() << stp1
<< "("<<tor1
<<")" << stp2
<< "("<<tor2
<<")";
1230 for (Int z
=0;z
<samples
;z
++) {
1232 x0
= randpi(ran
), xdot0
= randpi(ran
), y0
= randpi(ran
), ydot0
= randpi(ran
), phi0
= randpi(ran
);
1233 IC_setup(y
,x0
,xdot0
,y0
,ydot0
,phi0
,params
);
1234 bifparam(v
,param1
,stp1
,params
);
1235 bifparam(v
,param2
,stp2
,params
);
1236 // inner for: iterate over runge-kutta
1239 tau
= 2*PI
/params
[4];
1241 //qDebug() << "tau: " << tau << "\tL: " << L;
1243 xtest
= x
+h
; // test next x increment
1244 //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1245 if (xtest
>= k
*tau
) {
1246 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1247 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1249 //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1250 derivs(x
,y
,dydx
,params
);
1251 derivw(x
,y
, dW
,params
,gauss
);
1252 euler(y
,dydx
,dW
,x
,hnew
,yout
);
1255 derivs(x
,y
,dydx
,params
);
1256 derivw(x
,y
, dW
,params
,gauss
);
1257 euler(y
,dydx
,dW
,x
,h
,yout
); // xold
1258 for (Int j
=0;j
<n
;j
++)
1259 y
[j
] = yout
[j
]; // y[x+h] !!!!!!!
1260 x
+= h
; // do it for real
1262 Int m_x
, m_y
, m_phi
;
1263 m_x
= yout
[0]/L
; m_y
= yout
[1]/L
; m_phi
= yout
[2]/L
;
1264 avg
[0] = (Doub
) m_x
/k
; avg
[1] = (Doub
) m_y
/k
;
1265 avg
[2] = (Doub
) m_phi
/k
;
1266 set_color(avg
[0], vx_color
);
1267 set_color(avg
[1], vy_color
);
1268 out
<< stp1
<< "\t" << stp2
<< "\t"
1269 << avg
[0] << "\t" << vx_color
<< "\t" << avg
[1]
1270 << "\t" << vy_color
<< "\t" << avg
[2] << endl
;
1273 } // end 2. outer for
1274 } // end 1. outer for
1276 qDebug() << "moving" << name
<< "to safety";
1277 file
.copy(name
, file_rk_dualbif
);
1283 throw("No method specified!");
1286 rkValues_xdot
.clear();
1288 rkValues_ydot
.clear();
1289 rkValues_phi
.clear();
1290 rkValues_phidot
.clear();
1291 rkAvgvelo_x
.clear();
1292 rkAvgvelo_y
.clear();
1293 rkAvgvelo_phi
.clear();
1296 euValues_xdot
.clear();
1298 euValues_ydot
.clear();
1299 euValues_phi
.clear();
1300 euValues_phidot
.clear();
1301 euAvgvelo_x
.clear();
1302 euAvgvelo_y
.clear();
1303 euAvgvelo_phi
.clear();
1307 //=============================================================================