calculator.cpp

   1 //=============================================================================
   2 //
   3 // CLASS Calculator
   4 // contains all the magic math
   5 //
   6 //=============================================================================
   7
   8
   9 //== INCLUDES =================================================================
  10
  11 #include "calculator.h"
  12
  13
  14 //== IMPLEMENTATION ===========================================================
  15
  16
  17 Calculator::Calculator(QuPlot *parent) :
  18         QDialog(parent)
  19 {
  20         c = parent;
  21         derive_init();
  22
  23 }
  24
  25 Calculator::~Calculator()
  26 {
  27 }
  28
  29
  30 //-----------------------------------------------------------------------------
  31 void Calculator::derive_init()
  32 {
  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";
  39
  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";
  46
  47     file_rk_conv = "../rk_conv.dat";
  48     file_eu_conv = "../eu_conv.dat";
  49
  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";
  56
  57     file_rk_dualbif = "../rk_dualbif.dat";
  58     file_eu_dualbif = "../eu_dualbif.dat";
  59
  60     v = c->valuesMap.value(1);
  61
  62 }
  63
  64 void Calculator::derivs(const Doub x, VecDoub_I &y, VecDoub_O &dydx,
  65                         const Doub params[])
  66 {
  67     //params[0] = a;
  68     //params[1] = A;
  69     //params[2] = F;
  70     //params[3] = l;
  71     //params[4] = w;
  72     //params[5] = eta1;
  73     //params[6] = eta2;
  74     //params[7] = l1;
  75     //params[8] = l2;
  76     //params[9] = T;
  77     //params[10] = theta;
  78     //params[11] = e0;
  79     //params[12] = e1;
  80     //params[13] = e2;
  81     //params[14] = e3;
  82     //params[15] = fx;
  83     //params[16] = fy;
  84     //params[17] = psi;
  85     //params[18] = U;
  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];
  91
  92     switch (v) {
  93     case 1:
  94         G = sqrt(2.0*params[5]*params[9])*xi;
  95         dydx[0] = y[1];
  96         dydx[1] = -params[5]*y[1]-sin(y[0])+params[1]*sin(params[4]*x)+params[2]+G;
  97         break;
  98     case 2:
  99         dydx[0] = y[1];
 100         dydx[1] = -params[5]*y[1]-y[0];
 101         break;
 102     case 3:
 103         dydx[0] = y[1];
 104         dydx[1] = -y[0];
 105         break;
 106     case 4:
 107
 108         /*
 109         xix1 = gauss.dev();
 110         xix2 = gauss.dev();
 111         xiy1 = gauss.dev();
 112         xiy2 = gauss.dev();
 113         */
 114         /*
 115         eta = params[5]+params[6];
 116         e0 = 1.0/eta;
 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]));
 120
 121         l1 = params[3]*e0*params[5];
 122         l2 = params[3]*e0*params[6];
 123
 124         */
 125
 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];
 131
 132         //G = sqrt(2.0*eta*params[9]);
 133         G = 0;
 134
 135         c = cos(y[2]);
 136         s = sin(y[2]);
 137
 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;
 142
 143         st = sin(params[4]*x);
 144
 145         //qDebug() << st << "=" << "sin("<<params[4]<<"*"<<x<<")";
 146
 147         ax = cos(params[1])*params[0]*st;
 148         ay = sin(params[1])*params[0]*st;
 149
 150         //fx = cos(params[10])*params[2];
 151         //fy = sin(params[10])*params[2];
 152
 153         stx = ax+params[15];
 154         sty = ay+params[16];
 155
 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;
 162
 163
 164
 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;
 169
 170         cp = cos(params[17]);
 171         sp = sin(params[17]);
 172
 173         m1 = cp*x1-sp*y1;
 174         p1 = sp*x1+cp*y1;
 175         m2 = cp*x2-sp*y2;
 176         p2 = sp*x2+cp*y2;
 177
 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);
 182
 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))*/;
 186
 187         break;
 188     }
 189 }
 190
 191 void Calculator::derivw(const Doub x, VecDoub_I &y, VecDoub_O &dW,
 192        const Doub params[], Normaldev &gauss)
 193 {
 194     //params[0] = a;
 195     //params[1] = A;
 196     //params[2] = F;
 197     //params[3] = l;
 198     //params[4] = w;
 199     //params[5] = eta1;
 200     //params[6] = eta2;
 201     //params[7] = l1;
 202     //params[8] = l2;
 203     //params[9] = T;
 204     //params[10] = theta;
 205     //params[11] = e0;
 206     //params[12] = e1;
 207     //params[13] = e2;
 208     //params[14] = e3;
 209     //params[15] = fx;
 210     //params[16] = fy;
 211     //params[17] = psi;
 212     //params[18] = U;
 213
 214     Doub xix1, xix2, xiy1, xiy2;
 215     Doub c, s, T, l, eta1, eta2;
 216
 217     xix1 = gauss.dev();
 218     xix2 = gauss.dev();
 219     xiy1 = gauss.dev();
 220     xiy2 = gauss.dev();
 221
 222     T = params[9];
 223     l = params[7] + params[8];
 224     eta1 = params[5];
 225     eta2 = params[6];
 226
 227     c = cos(y[2]);
 228     s = sin(y[2]);
 229
 230     switch (v) {
 231         case 4:
 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));
 235         break;
 236     }
 237 }
 238 //-----------------------------------------------------------------------------
 239
 240
 241 void Calculator::rk4(VecDoub_I &y, VecDoub_I &dydx, const Doub x, const Doub h,
 242                      VecDoub_O &yout, const Doub params[])
 243 {
 244     Int n = y.size();
 245     VecDoub dym(n),dyt(n),yt(n);
 246     Doub hh = h*0.5;
 247     Doub h6 = h/6.0;
 248     Doub xh = x+hh;
 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];
 255         dym[i] += dyt[i];
 256     }
 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]);
 260     }
 261 }
 262
 263
 264 void Calculator::euler(VecDoub_I &y, VecDoub_I &dydx, VecDoub_I &dW, const Doub x,
 265                        const Doub h, VecDoub_O &yout)
 266 {
 267     Int n = y.size();
 268     for (Int i=0;i<n;i++) {
 269         yout[i] = y[i]+h*dydx[i]+sqrt(h)*dW[i];
 270     }
 271 }
 272 //-----------------------------------------------------------------------------
 273
 274
 275 /*
 276  * returns a random number between 0 and 2*pi
 277  */
 278 double Calculator::randpi(Ran &r)
 279 {
 280     Doub s = 2.0*M_PI*r.doub();
 281     return s;
 282 }
 283
 284 //-----------------------------------------------------------------------------
 285
 286 /*
 287  * Writes the of two maps into a tab-separated file
 288  * The second map *must* be specified.
 289  */
 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,
 294                                const Doub params[])
 295 {
 296     qDebug() << "writing to file" << name;
 297
 298
 299     QFile file(name);
 300     QTextStream out(&file);
 301     file.open(QIODevice::WriteOnly);
 302     out << "#";
 303     for (Int i=0;i<10;i++)
 304         out << "params["<<i<<"] = " <<  params[i] << " ";
 305     out << endl;
 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);
 310     switch(v) {
 311         case 1:
 312         while (x.hasNext()) {
 313             x.next();
 314             out << x.key() << "\t" << x.value() << endl;
 315         }
 316         break;
 317         case 4:
 318         if (map_x.size() > 0 && map_y.size() == 0 && map_phi.size() == 0) {
 319             while (x.hasNext()) {
 320                 x.next();
 321                 out << x.key() << "\t" << x.value() << endl;
 322             }
 323         }
 324         else if (map_x.size() > 0 && map_y.size() > 0 && map_phi.size() > 0) {
 325             while (x.hasNext() && y.hasNext() && phi.hasNext()) {
 326                 x.next();
 327                 y.next();
 328                 phi.next();
 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;
 331             }
 332         }
 333         break;
 334     }
 335     file.close();
 336     qDebug () << "done.";
 337 }
 338
 339
 340 //-----------------------------------------------------------------------------
 341
 342
 343 void Calculator::diffusion(VecDoub_I &msd, const Int & samples, Doub & sigma)
 344 {
 345     Doub mu,sum;
 346     VecDoub var(samples,0.0);
 347
 348     sum = 0;
 349     for(Int i=0;i<samples;i++) {
 350         //qDebug() << "msd["<<i<<"]" << msd[i];
 351         sum += msd[i];
 352     }
 353     mu = sum/samples;
 354     //qDebug() << "mu:" << mu;
 355
 356     for(Int i=0;i<samples;i++) {
 357         var[i] = SQR(msd[i])-SQR(mu);
 358         //qDebug() << "var["<<i<<"]" << var[i];
 359     }
 360     sum = 0;
 361     for(Int j=0;j<samples;j++) {
 362         sum += var[j];
 363     }
 364
 365     sigma = sum/samples;
 366
 367     //qDebug() << "sigma:" << sigma;
 368 }
 369
 370 /*
 371  * This function is kinda redundant right now, but will be needed later
 372  */
 373 void Calculator::gen_gnuplot(const QString & name)
 374 {
 375
 376     qDebug() << "generating gnuplot file";
 377
 378     QFile file(name);
 379     QTextStream out(&file);
 380     file.open(QIODevice::WriteOnly);
 381
 382     out << "L = 2*pi;";
 383     out << "unset key;";
 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;";
 392     out << "set xtics;";
 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;";
 413
 414 }
 415
 416 void Calculator::filter(QMap<double, double> & map, const Doub & eps)
 417 {
 418     qDebug() << "filter crap smaller than " << eps;
 419     Int size, k;
 420     Doub diff, key;
 421     QList<double> values, list;
 422     list = map.uniqueKeys();
 423     //qDebug() << list;
 424     size = list.size();
 425     //qDebug() << "size of map:" << size;
 426     for(Int i = 0; i < size; i++) {
 427         key = list.at(i);
 428         //qDebug() << "["<<key<<"]" << "writing values into list";
 429         values = map.values(key);
 430         //qDebug() << "["<<key<<"]" << values.size() << "values for" << key;
 431         k = 0;
 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;
 435             if (diff <= eps) {
 436                 //qDebug() << diff << "<=" << eps;
 437                 //qDebug() << "["<<k<<"]" << "removing" << values.at(k) << "from list";
 438                 values.replace(k+1,values.at(k));
 439             }
 440         } // end for
 441         for (int j = 0; j < values.size(); j++)
 442             //qDebug() << "["<<j<<"]" << values.at(j);
 443         //qDebug() << "removing" << key << "from map";
 444         map.remove(key);
 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));
 448         }
 449         //qDebug() << map.values(key);
 450         //qDebug() << endl << "end while" << endl;
 451     }
 452 }
 453
 454 void  Calculator::IC_setup(VecDoub &y, Doub x0, Doub y0, Doub xdot0, Doub ydot0,
 455                           Doub phi0, Doub params[])
 456 {
 457     //Int n = y.size();
 458     const Int sys = c->valuesMap.value(1);
 459
 460     /*
 461     qDebug() << "using " << "\tx0 = " << x0 << "\ty0 = " << y0
 462             << "\txdot0 = " << xdot0 << "\tydot0 = " << ydot0
 463             << "\tphi0 = " << fmod(phi0,360) << endl;
 464     */
 465
 466     switch(sys) {
 467         case 1:
 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;
 474         params[6] = 0;
 475         params[9] = c->valuesMap.value(32); // T;
 476         params[10] = 0;
 477         y[0] = x0, y[1] = xdot0;
 478         break;
 479         case 4:
 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
 498         params[12] = e1;
 499         params[13] = e2;
 500         params[14] = e3;
 501         params[15] = fx;
 502         params[16] = fy;
 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;
 506         break;
 507     }
 508     /*
 509     for (Int i = 0; i<n;i++)
 510         qDebug() << "y("<<i<<") =" << y[i];
 511     qDebug() << "eta1 = " << eta1 << "\teta2 = " << eta2 << "\tA = " << A << "\tF = " << F;
 512     */
 513 }
 514
 515 void Calculator::bifparam(const int &v, const QString &param, const Doub &r,
 516                           Doub params[])
 517 {
 518     switch(v) {
 519         case 1:
 520         if (param == "E")
 521             params[5] = r;
 522         else if (param == "F")
 523             params[2] = r;
 524         else if (param == "A")
 525             params[1] = r;
 526         else
 527             throw ("Fatal error!");
 528         break;
 529         case 4:
 530         if (param == "A")
 531             params[1] = fmod(r,360)*PI/180;
 532         else if (param == "F") {
 533             params[2] = r;
 534             params[15] = cos(params[10])*params[2];
 535             params[16] = sin(params[10])*params[2];
 536         }
 537         else if (param == "a")
 538             params[0] = r;
 539         else if (param == "w")
 540             params[4] = r;
 541         else if (param == "l") {
 542             params[3] = r;
 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;
 549         }
 550         else if (param == "e") {
 551             params[5] = r;
 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;
 558         }
 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];
 563         }
 564         else
 565             throw("Fatal error!");
 566         break;
 567         qDebug() << "Bifurcation parameter is " << param;
 568     }
 569 }
 570
 571 void Calculator::set_color(const Doub vel, QString & color)
 572 {
 573     Doub ctol = 0.05;
 574     if (vel > 1.0-ctol && vel < 1.0+ctol)
 575         color = "red";
 576     else if (vel > 0.5-ctol && vel < 0.5+ctol)
 577         color = "orange";
 578     else if (vel > 0.0-ctol && vel < 0.0+ctol)
 579         color = "grey";
 580     else if (vel > -0.5-ctol && vel < -0.5+ctol)
 581         color = "green";
 582     else if (vel > -1.0-ctol && vel < -1.0+ctol)
 583         color = "blue";
 584     else color = "black";
 585 }
 586
 587 void Calculator::calc2(const Int n)
 588 {
 589     // setting up a few start values
 590     qDebug() << "calc2 has been called: " << "calc2("<<n<<")";
 591     Normaldev gauss(0,1,time(NULL));
 592     Ran ran(time(NULL));
 593     Doub x1,x2,y1,y2,xtest;
 594     Doub L, tau, sigma;
 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);
 618     Doub params[256];
 619     //const Doub L = 1;
 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;
 625
 626     Doub x, hnew;
 627     Doub trajectory_length, deflection, velocity_angle, bias_angle;
 628     Int k, hit;
 629
 630     Doub xpi[samples][s];
 631     Doub ypi[samples][s];
 632     Doub phipi[samples][s];
 633     VecDoub y(n,0.0);
 634     VecDoub dydx(n,0.0);
 635     VecDoub dW(n, 0.0);
 636     VecDoub yout(n,0.0);
 637     VecDoub avg(n,0.0);
 638
 639     if (ruk) {
 640         switch(M) {
 641             case 0:
 642             qDebug() << "standard mode";
 643             qDebug() << "stepsize:" << h;
 644             qDebug() << "number of steps:" << s;
 645             x = t0;
 646             IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 647             k = 0;
 648             L = 2*PI;
 649             tau = 2*PI/params[4];
 650             while(k < s) {
 651                 xtest = x+h;
 652                 if (xtest >= k*tau) {
 653                     k++;
 654                 }
 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++)
 658                     y[j] = yout[j];
 659                 x += h;
 660                 if (k >= s*fromt) {
 661                     switch(v) {
 662                         case 1:
 663                         rkValues_x.insert(x,yout[0]);
 664                         break;
 665                         case 4:
 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);
 681                         break;
 682                     }
 683                 }
 684             } // end while
 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);
 691             break;
 692             case 1:
 693             qDebug() << "converge mode";
 694             Int iter;
 695             Doub hh, ss;
 696             for (iter=1;from+(iter*hstep)<=to;iter++) {
 697                 for (Int z=0;z<samples;z++) {
 698                     qDebug() << "at step" << z << "of" << samples;
 699                     x = t0;
 700                     hh = from+(iter*hstep);
 701                     ss = (s/hh);
 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);
 705                     k = 0;
 706                     L = 2*PI;
 707                     tau = 2*PI/params[4];
 708                     while(k < s) {
 709                         //qDebug() << k << "<" << s;
 710                         //qDebug() << k << tau << k*tau;
 711                         xtest = x+hh;
 712                         if (xtest >= k*tau) {
 713                             //qDebug() << xtest << ">=" << k << "*" << tau << "("<<k*tau<<")";
 714                             hnew = k*tau-x;
 715                             derivs(x,y,dydx,params);
 716                             rk4(y,dydx,x,hnew,yout,params);
 717                             k++;
 718                             x += hnew;
 719                             //qDebug() << x;
 720                         } else {
 721                             //qDebug() << xtest << "=" << x<<"+"<<hh;
 722                             derivs(x,y,dydx,params);
 723                             rk4(y,dydx,x,hh,yout,params); // xold
 724                             x += hh;
 725                         }
 726                         for (Int j=0;j<n;j++)
 727                             y[j] = yout[j];
 728                     } // end while
 729                     switch(v) {
 730                         case 1:
 731                         rkValues_x.insert(hh,yout[0]);
 732                         break;
 733                         case 4:
 734                         rkValues_x.insert(hh,yout[0]);
 735                         rkValues_y.insert(hh,yout[1]);
 736                         rkValues_phi.insert(hh,yout[2]);
 737                         break;
 738                     }
 739                     x = t0;
 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);
 742                     k = 0;
 743                     L = 2*PI;
 744                     tau = 2*PI/params[4];
 745                     while (k < s) {
 746                         xtest = x+hh;
 747                         if (xtest >= k*tau) {
 748                             hnew = k*tau-x;
 749                             derivs(x,y,dydx,params);
 750                             derivw(x,y,dW,params,gauss);
 751                             euler(y,dydx,dW,x,hnew,yout);
 752                             k++;
 753                             x += hnew;
 754                         } else {
 755                             derivs(x,y,dydx,params);
 756                             derivw(x,y,dW,params,gauss);
 757                             euler(y,dydx,dW,x,hh,yout);
 758                             x += hh;
 759                         }
 760                         for (Int j=0;j<n;j++)
 761                             y[j] = yout[j];
 762                     } // end while
 763                     switch(v) {
 764                         case 1:
 765                         euValues_x.insert(hh,yout[0]);
 766                         break;
 767                         case 4:
 768                         euValues_x.insert(hh,yout[0]);
 769                         euValues_y.insert(hh,yout[1]);
 770                         euValues_phi.insert(hh,yout[2]);
 771                         break;
 772                     } // end switch
 773                 } // end samples
 774             } // end for
 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);
 778             break;
 779             case 2:
 780             qDebug() << "bifurcation mode";
 781             qDebug() << "bifurcation parameter: " << param1;
 782             qDebug() << "Number of samples : " << samples;
 783             Doub r;
 784             for (Int j=0;(j*stepr1)<=(tor1-fromr1);j++) {
 785                 r = fromr1+j*stepr1;
 786                 qDebug() << "at step" << j << "of" << (tor1-fromr1)/stepr1;
 787                 for (Int z=0;z<samples;z++) {
 788                     x = t0;
 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
 797                     k = 0;
 798                     L = 2*PI;
 799                     tau = 2*PI/params[4];
 800                     Doub hnew;
 801                     //qDebug() << "tau: " << tau << "\tL: " << L;
 802                     while (k<s) {
 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;
 808                             hnew = k*tau-x;
 809                             //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 810                             derivs(x,y,dydx,params);
 811                             rk4(y,dydx,x,hnew,yout,params);
 812                             switch(v) {
 813                                 case 1:
 814                                 xpi[z][k] = ABS(fmod(yout[0],(2*PI)));
 815                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 816                                 if (k >= s*fromt) {
 817                                     rkValues_x.insert(r,xpi[z][k]);
 818                                     rkValues_xdot.insert(r,dydx[0]);
 819                                     }
 820                                 break;
 821                                 case 4:
 822                                 xpi[z][k] = yout[0];
 823                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 824                                 ypi[z][k] = yout[1];
 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));
 829                                 if (k >= s*fromt) {
 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]);
 836                                 }
 837                                 break;
 838                             } // end switch
 839                             k++;
 840                         }
 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
 846                     } // end while
 847                     //qDebug() << "at" << x << "in time for k = " << k;
 848                     switch(v) {
 849                         Int m_x, m_y, m_phi;
 850                         case 1:
 851                         m_x = yout[0]/L;
 852                         avg[0] = (Doub) m_x/k;
 853                         rkAvgvelo_x.insert(r,avg[0]);
 854                         break;
 855                         case 4:
 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) {
 861                             velocity_angle = PI;
 862                         } else {
 863                             velocity_angle = (acos((yout[0]-x0)/trajectory_length));
 864                         }
 865                         bias_angle = params[10];
 866                         deflection = (180/PI)*ABS(velocity_angle - bias_angle);
 867                         /*
 868                         qDebug() << velocity_angle << "("<<yout[0]-x0<<")"
 869                                 << "-" << bias_angle << "="
 870                                 << deflection
 871                                 << "[ length:"<<trajectory_length<<"]";
 872                         */
 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);
 877                         break;
 878                     } // end switch
 879                 } // samples loop
 880                 VecDoub msd_x(samples,0.0);
 881                 VecDoub msd_y(samples,0.0);
 882                 VecDoub msd_phi(samples,0.0);
 883                 Doub sigma_x[s-1];
 884                 Doub sigma_y[s-1];
 885                 Doub sigma_phi[s-1];
 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];
 895                     }
 896                     diffusion(msd_x, samples, sigma);
 897                     sigma_x[i] = sigma;
 898                     diffusion(msd_y, samples, sigma);
 899                     sigma_y[i] = sigma;
 900                     diffusion(msd_phi, samples, sigma);
 901                     sigma_phi[i] = sigma;
 902                     //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
 903                 }
 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]);
 909                 }
 910                //qDebug() << rkVariance_x;
 911             } // end outer for
 912             /*
 913              * we need to do something about this
 914              */
 915             /*
 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);
 933             break;
 934             case 3:
 935             qDebug() << "dual bifurcation mode";
 936             qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
 937             qDebug() << "Number of samples : " << samples;
 938             Doub stp1, stp2;
 939             QString vx_color, vy_color;
 940             QString name = "tmp.dat";
 941             QFile file(name);
 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++) {
 951                         x = t0;
 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);
 956                         /*
 957                         out << "# ";
 958                         for (Int i=0;i<10;i++)
 959                             out << "params["<<i<<"] = " <<  params[i] << " ";
 960                         out << endl;
 961                         */
 962                         // inner for: iterate over runge-kutta
 963                         k = 0;
 964                         L = 2*PI;
 965                         tau = 2*PI/params[4];
 966                         Doub hnew;
 967                         //qDebug() << "tau: " << tau << "\tL: " << L;
 968                         while (k<s) {
 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;
 974                                 hnew = k*tau-x;
 975                                 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 976                                 derivs(x,y,dydx,params);
 977                                 rk4(y,dydx,x,hnew,yout,params);
 978                                 k++;
 979                             }
 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
 985                         } // end while
 986                         Int m_x, m_y, m_phi;
 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;
 995                     } // end samples
 996                 } // end 2. outer for
 997             } // end 1. outer for
 998             file.close();
 999             qDebug() << "moving" << name << "to safety";
1000             file.copy(name, file_rk_dualbif);
1001             break;
1002         } // end switch
1003     } // end if
1004
1005     if (eul) {
1006         switch(M) {
1007           case 0:
1008           qDebug() << "(euler) standard mode";
1009           qDebug() << "stepsize:" << h;
1010           qDebug() << "number of steps:" << s;
1011           x = t0;
1012           IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1013           k = 0;
1014           hit = 0;
1015           L = 2*PI;
1016           tau = 2*PI/params[4];
1017           while(k < s) {
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];
1027               if (x >= hit*0.1) {
1028                   //qDebug() << x <<"=="<< hit*0.1;
1029                   //qDebug() << "writing2 to map";
1030                   hit++;
1031                   //qDebug() << "next value to write out" << hit*0.1;
1032                   if (k >= s*fromt) {
1033                       switch(v) {
1034                       case 1:
1035                           euValues_x.insert(x,yout[0]);
1036                       break;
1037                       case 4:
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);
1053                         break;
1054                     } // end fromt switch
1055                   } // end fromt if
1056               } // end hit if
1057               //qDebug() << "making step of size" << h;
1058               for (Int j=0;j<n;j++)
1059                   y[j] = yout[j];
1060               x += h;
1061               //qDebug() << "new x is now" << x;
1062               if (x >= k*tau) {
1063                   //qDebug() << x << ">=" << k*tau;
1064                   //qDebug() << "incrementing k";
1065                   k++;
1066                   //qDebug() << "new value of x needed is" << k*tau;
1067               }
1068           } // end while
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);
1075           break;
1076           case 2:
1077           qDebug() << "bifurcation mode";
1078           qDebug() << "bifurcation parameter: " << param1;
1079           qDebug() << "Number of samples : " << samples;
1080           Doub r;
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++) {
1085                   x = t0;
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
1091                   k = 0;
1092                   L = 2*PI;
1093                   tau = 2*PI/params[4];
1094                   Doub hnew;
1095                   //qDebug() << "tau: " << tau << "\tL: " << L;
1096                   while (k<s) {
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;
1102                           hnew = k*tau-x;
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);
1107                           switch(v) {
1108                           case 1:
1109                               xpi[z][k] = yout[0];
1110                               if (k >= s*fromt) {
1111                                   euValues_x.insert(r,ABS(fmod(xpi[z][k],L)));
1112                                   euValues_xdot.insert(r,dydx[0]);
1113                               }
1114                               break;
1115                           case 4:
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));
1123                               if (k >= s*fromt) {
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]);
1130                               }
1131                               break;
1132                           } // end switch
1133                           k++;
1134                       }
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
1141                   } // end while
1142                   //qDebug() << "at" << x << "in time for k = " << k;
1143                   switch(v) {
1144                       Int m_x, m_y, m_phi;
1145                       case 1:
1146                       m_x = yout[0]/L;
1147                       avg[0] = (Doub) m_x/k;
1148                       euAvgvelo_x.insert(r,avg[0]);
1149                       break;
1150                       case 4:
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;
1156                       } else {
1157                           velocity_angle = (acos((yout[0]-x0)/trajectory_length));
1158                       }
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);
1165                       break;
1166                   } // end switch
1167               } // samples loop
1168               VecDoub msd_x(samples,0.0);
1169               VecDoub msd_y(samples,0.0);
1170               VecDoub msd_phi(samples,0.0);
1171               Doub sigma_x[s-1];
1172               Doub sigma_y[s-1];
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];
1183                   }
1184                   diffusion(msd_x, samples, sigma);
1185                   sigma_x[i] = sigma;
1186                   diffusion(msd_y, samples, sigma);
1187                   sigma_y[i] = sigma;
1188                   diffusion(msd_phi, samples, sigma);
1189                   sigma_phi[i] = sigma;
1190                   //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
1191               }
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]);
1197               }
1198               //qDebug() << euVariance_x;
1199           } // end outer for
1200           /*
1201            * we need to do something about this
1202            */
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);
1213           break;
1214           case 3:
1215           qDebug() << "dual bifurcation mode";
1216           qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
1217           qDebug() << "Number of samples : " << samples;
1218           Doub stp1,stp2;
1219           QString vx_color, vy_color;
1220           QString name = "tmp.dat";
1221           QFile file(name);
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++) {
1231                       x = t0;
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
1237                       Int k = 0;
1238                       L = 2*PI;
1239                       tau = 2*PI/params[4];
1240                       Doub hnew;
1241                       //qDebug() << "tau: " << tau << "\tL: " << L;
1242                       while (k<s) {
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;
1248                               hnew = k*tau-x;
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);
1253                               k++;
1254                           }
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
1261                       } // end while
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;
1271
1272                   } // end samples
1273               } // end 2. outer for
1274           } // end 1. outer for
1275           file.close();
1276           qDebug() << "moving" << name << "to safety";
1277           file.copy(name, file_rk_dualbif);
1278           break;
1279         } // end switch
1280     } // end if
1281
1282     if (!ruk && !eul)
1283         throw("No method specified!");
1284
1285     rkValues_x.clear();
1286     rkValues_xdot.clear();
1287     rkValues_y.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();
1294
1295     euValues_x.clear();
1296     euValues_xdot.clear();
1297     euValues_y.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();
1304
1305 } // end calculator
1306
1307 //=============================================================================