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     file_rk_vel = "../rk_velo.dat";
  40
  41     file_eu_bif = "../eu_bifurcate.dat";
  42     file_eu_avg = "../eu_bif_avg_velo.dat";
  43     file_eu_out = "../eu_out.dat";
  44     file_eu_pt1 = "../eu_point1.dat";
  45     file_eu_pt2 = "../eu_point2.dat";
  46     file_eu_biv = "../eu_bif_velo.dat";
  47     file_eu_vel = "../eu_velo.dat";
  48
  49     file_rk_conv = "../rk_conv.dat";
  50     file_eu_conv = "../eu_conv.dat";
  51
  52     file_rk_variance = "../rk_variance.dat";
  53     file_rk_distance = "../rk_msd.dat";
  54     file_rk_defangle = "../rk_defangle.dat";
  55     file_eu_variance = "../eu_variance.dat";
  56     file_eu_distance = "../eu_msd.dat";
  57     file_eu_defangle = "../eu_defangle.dat";
  58
  59     file_rk_dualbif = "../rk_dualbif.dat";
  60     file_eu_dualbif = "../eu_dualbif.dat";
  61
  62     v = c->valuesMap.value(1);
  63
  64 }
  65
  66 void Calculator::derivs(const Doub x, VecDoub_I &y, VecDoub_O &dydx,
  67                         const Doub params[])
  68 {
  69     //params[0] = a;
  70     //params[1] = A;
  71     //params[2] = F;
  72     //params[3] = l;
  73     //params[4] = w;
  74     //params[5] = eta1;
  75     //params[6] = eta2;
  76     //params[7] = l1;
  77     //params[8] = l2;
  78     //params[9] = T;
  79     //params[10] = theta;
  80     //params[11] = e0;
  81     //params[12] = e1;
  82     //params[13] = e2;
  83     //params[14] = e3;
  84     //params[15] = fx;
  85     //params[16] = fy;
  86     //params[17] = psi;
  87     //params[18] = U;
  88     //static Doub G, xi, xi1, xi2, c, s, st, ax, ay, e0, e1, e2, e3, eta, l1, l2;
  89     static Doub xi, G, s, c, st, ax, ay, cp, sp, m1, m2, p1, p2;
  90     static Doub x1, y1, x2, y2, F1x, F1y, F2x, F2y, stx, sty;
  91     //for (Int i = 0;i<10;i++)
  92      //   qDebug() << "params["<<i<<"] :" << params[i];
  93
  94     switch (v) {
  95     case 1:
  96         G = sqrt(2.0*params[5]*params[9])*xi;
  97         dydx[0] = y[1];
  98         dydx[1] = -params[5]*y[1]-sin(y[0])+params[1]*sin(params[4]*x)+params[2]+G;
  99         break;
 100     case 2:
 101         dydx[0] = y[1];
 102         dydx[1] = -params[5]*y[1]-y[0];
 103         break;
 104     case 3:
 105         dydx[0] = y[1];
 106         dydx[1] = -y[0];
 107         break;
 108     case 4:
 109
 110         /*
 111         xix1 = gauss.dev();
 112         xix2 = gauss.dev();
 113         xiy1 = gauss.dev();
 114         xiy2 = gauss.dev();
 115         */
 116         /*
 117         eta = params[5]+params[6];
 118         e0 = 1.0/eta;
 119         e1 = 1.0/(params[5]*params[3]);
 120         e2 = 1.0/(params[6]*params[3]);
 121         e3 = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 122
 123         l1 = params[3]*e0*params[5];
 124         l2 = params[3]*e0*params[6];
 125
 126         */
 127
 128        //qDebug() << eta << "=" << params[5] << "+" << params[6];
 129        //qDebug() << e0 << "=" << "1.0/"<< eta;
 130        //qDebug() << e1 << "=" << "1.0/"<<params[5]<<"*"<<params[3];
 131        //qDebug() << e2 << "=" << "1.0/"<<params[6]<<"*"<<params[3];
 132        //qDebug() << e3 << "=" << "1.0/"<<params[5]<<"-"<<"1.0/"<<params[6];
 133
 134         //G = sqrt(2.0*eta*params[9]);
 135         G = 0;
 136
 137         c = cos(y[2]);
 138         s = sin(y[2]);
 139
 140         x1 = y[0] + params[7]*c;
 141         y1 = y[1] + params[7]*s;
 142         x2 = y[0] - params[8]*c;
 143         y2 = y[1] - params[8]*s;
 144
 145         st = sin(params[4]*x);
 146
 147         //qDebug() << st << "=" << "sin("<<params[4]<<"*"<<x<<")";
 148
 149         ax = cos(params[1])*params[0]*st;
 150         ay = sin(params[1])*params[0]*st;
 151
 152         //fx = cos(params[10])*params[2];
 153         //fy = sin(params[10])*params[2];
 154
 155         stx = ax+params[15];
 156         sty = ay+params[16];
 157
 158         //qDebug() << "ax =" << ax << "=" << params[0] << "*" << "cos("<<params[1]<<")";
 159         //qDebug() << "ay =" << ay << "=" << params[0] << "*" << "sin("<<params[1]<<")";
 160         //qDebug() << "fx =" << fx << "=" << "cos("<<params[10]<<")*"<<params[2];
 161         //qDebug() << "fy =" << fy << "=" << "sin("<<params[10]<<")*"<<params[2];
 162         //qDebug() << "stx =" << stx <<"="<< ax << "*" << st << "+"<< fx;
 163         //qDebug() << "sty =" << sty <<"="<< ay << "*" << st << "+" << fy;
 164
 165
 166
 167         //qDebug() << y[0] + params[7]*c << "=" << y[0] << "+" << params[7] << "*" << c;
 168         //qDebug() << y[1] + params[7]*c << "=" << y[1] << "+" << params[7] << "*" << s;
 169         //qDebug() << y[0] - params[8]*c << "=" << y[0] << "-" << params[8] << "*" << c;
 170         //qDebug() << y[1] - params[8]*c << "=" << y[1] << "-" << params[8] << "*" << s;
 171
 172         cp = cos(params[17]);
 173         sp = sin(params[17]);
 174
 175         m1 = cp*x1-sp*y1;
 176         p1 = sp*x1+cp*y1;
 177         m2 = cp*x2-sp*y2;
 178         p2 = sp*x2+cp*y2;
 179
 180         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);
 181         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);
 182         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);
 183         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);
 184
 185         dydx[0] = (F1x+F2x+2.0*stx/*+sqrt(2*params[9])*(xix1+xix2)*/)*params[11];
 186         dydx[1] = (F1y+F2y+2.0*sty/*+sqrt(2*params[9])*(xiy1+xiy2)*/)*params[11];
 187         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))*/;
 188
 189         break;
 190     }
 191 }
 192
 193 void Calculator::derivw(const Doub x, VecDoub_I &y, VecDoub_O &dW,
 194        const Doub params[], Normaldev &gauss)
 195 {
 196     //params[0] = a;
 197     //params[1] = A;
 198     //params[2] = F;
 199     //params[3] = l;
 200     //params[4] = w;
 201     //params[5] = eta1;
 202     //params[6] = eta2;
 203     //params[7] = l1;
 204     //params[8] = l2;
 205     //params[9] = T;
 206     //params[10] = theta;
 207     //params[11] = e0;
 208     //params[12] = e1;
 209     //params[13] = e2;
 210     //params[14] = e3;
 211     //params[15] = fx;
 212     //params[16] = fy;
 213     //params[17] = psi;
 214     //params[18] = U;
 215
 216     Doub xix1, xix2, xiy1, xiy2;
 217     Doub c, s, T, l, eta1, eta2;
 218
 219     xix1 = gauss.dev();
 220     xix2 = gauss.dev();
 221     xiy1 = gauss.dev();
 222     xiy2 = gauss.dev();
 223
 224     T = params[9];
 225     l = params[7] + params[8];
 226     eta1 = params[5];
 227     eta2 = params[6];
 228
 229     c = cos(y[2]);
 230     s = sin(y[2]);
 231
 232     switch (v) {
 233         case 4:
 234         dW[0] = sqrt(2*T)*(sqrt(eta1)*xix1+sqrt(eta2)*xix2);
 235         dW[1] = sqrt(2*T)*(sqrt(eta1)*xiy1+sqrt(eta2)*xiy2);
 236         dW[2] = (sqrt(2*T)/l)*((1/sqrt(eta1))*(-s*xix1+c*xiy1)-(1/sqrt(eta2))*(-s*xix2+c*xiy2));
 237         break;
 238     }
 239 }
 240 //-----------------------------------------------------------------------------
 241
 242
 243 void Calculator::rk4(VecDoub_I &y, VecDoub_I &dydx, const Doub x, const Doub h,
 244                      VecDoub_O &yout, const Doub params[])
 245 {
 246     Int n = y.size();
 247     VecDoub dym(n),dyt(n),yt(n);
 248     Doub hh = h*0.5;
 249     Doub h6 = h/6.0;
 250     Doub xh = x+hh;
 251     for (Int i=0;i<n;i++) yt[i] = y[i]+hh*dydx[i];
 252     derivs(xh,yt,dyt,params);
 253     for (Int i=0;i<n;i++) yt[i] = y[i]+hh*dyt[i];
 254     derivs(xh,yt,dym,params);
 255     for (Int i=0;i<n;i++) {
 256         yt[i] = y[i]+h*dym[i];
 257         dym[i] += dyt[i];
 258     }
 259     derivs(x+h,yt,dyt,params);
 260     for (Int i=0;i<n;i++) {
 261         yout[i] = y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
 262     }
 263 }
 264
 265
 266 void Calculator::euler(VecDoub_I &y, VecDoub_I &dydx, VecDoub_I &dW, const Doub x,
 267                        const Doub h, VecDoub_O &yout)
 268 {
 269     Int n = y.size();
 270     for (Int i=0;i<n;i++) {
 271         yout[i] = y[i]+h*dydx[i]+sqrt(h)*dW[i];
 272     }
 273 }
 274 //-----------------------------------------------------------------------------
 275
 276
 277 /*
 278  * returns a random number between 0 and 2*pi
 279  */
 280 double Calculator::randpi(Ran &r)
 281 {
 282     Doub s = 2.0*M_PI*r.doub();
 283     return s;
 284 }
 285
 286 //-----------------------------------------------------------------------------
 287
 288 /*
 289  * Writes the of two maps into a tab-separated file
 290  * The second map *must* be specified.
 291  */
 292 void Calculator::write_to_file(const QString & name,
 293                                const QMap<double, double> & map_x,
 294                                const QMap<double, double> & map_y,
 295                                const QMap<double, double> & map_phi,
 296                                const Doub params[])
 297 {
 298     qDebug() << "writing to file" << name;
 299
 300
 301     QFile file(name);
 302     QTextStream out(&file);
 303     file.open(QIODevice::WriteOnly);
 304     out << "#";
 305     for (Int i=0;i<10;i++)
 306         out << "params["<<i<<"] = " <<  params[i] << " ";
 307     out << endl;
 308     out << "r\tx\ty\tphi" << endl;
 309     QMapIterator<double,double> x(map_x);
 310     QMapIterator<double,double> y(map_y);
 311     QMapIterator<double,double> phi(map_phi);
 312     switch(v) {
 313         case 1:
 314         while (x.hasNext()) {
 315             x.next();
 316             out << x.key() << "\t" << x.value() << endl;
 317         }
 318         break;
 319         case 4:
 320         if (map_x.size() > 0 && map_y.size() == 0 && map_phi.size() == 0) {
 321             while (x.hasNext()) {
 322                 x.next();
 323                 out << x.key() << "\t" << x.value() << endl;
 324             }
 325         }
 326         else if (map_x.size() > 0 && map_y.size() > 0 && map_phi.size() > 0) {
 327             while (x.hasNext() && y.hasNext() && phi.hasNext()) {
 328                 x.next();
 329                 y.next();
 330                 phi.next();
 331                 //qDebug() << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value();
 332                 out << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value() << endl;
 333             }
 334         }
 335         break;
 336     }
 337     file.close();
 338     qDebug () << "done.";
 339 }
 340
 341
 342 //-----------------------------------------------------------------------------
 343
 344
 345 void Calculator::diffusion(VecDoub_I &msd, const Int & samples, Doub & sigma)
 346 {
 347     Doub mu,sum;
 348     VecDoub var(samples,0.0);
 349
 350     sum = 0;
 351     for(Int i=0;i<samples;i++) {
 352         //qDebug() << "msd["<<i<<"]" << msd[i];
 353         sum += msd[i];
 354     }
 355     mu = sum/samples;
 356     //qDebug() << "mu:" << mu;
 357
 358     for(Int i=0;i<samples;i++) {
 359         var[i] = SQR(msd[i])-SQR(mu);
 360         //qDebug() << "var["<<i<<"]" << var[i];
 361     }
 362     sum = 0;
 363     for(Int j=0;j<samples;j++) {
 364         sum += var[j];
 365     }
 366
 367     sigma = sum/samples;
 368
 369     //qDebug() << "sigma:" << sigma;
 370 }
 371
 372 /*
 373  * This function is kinda redundant right now, but will be needed later
 374  */
 375 void Calculator::gen_gnuplot(const QString & name)
 376 {
 377
 378     qDebug() << "generating gnuplot file";
 379
 380     QFile file(name);
 381     QTextStream out(&file);
 382     file.open(QIODevice::WriteOnly);
 383
 384     out << "L = 2*pi;";
 385     out << "unset key;";
 386     out << "set lmargin 10;";
 387     out << "set rmargin 3;";
 388     out << "set multiplot;";
 389     out << "set tmargin = 3;";
 390     out << "set bmargin = 0;";
 391     out << "set border 2+4+8;";
 392     out << "set size 1,0.66667;";
 393     out << "set origin 0.0,0.33333;";
 394     out << "set xtics;";
 395     out << "set format x \"\";";
 396     out << "unset xlabel;";
 397     out << "set ylabel \"x(kT) mod L \";";
 398     out << "set ytics (\"0\" 0, \"L/2\", \"L\" L);";
 399     out << "set yrange [0:L];";
 400     out << "plot \""<<file_rk_bif<<"\" w dots 1;";
 401     out << "unset title;";
 402     out << "set tmargin 0;";
 403     out << "set bmargin 3;";
 404     out << "unset format;";
 405     out << "set size 1,0.33333;";
 406     out << "set border 1+2+4+8;";
 407     out << "set origin 0.0,0.0;";
 408     out << "set xtics 0,0.1;";
 409     out << "set xlabel \"a\";";
 410     out << "set ylabel \"v\";";
 411     out << "set yrange [-0.75:0.75];";
 412     out << "set ytics -0.5,0.5;";
 413     out << "plot \""<<file_rk_avg<<"\" w dots 3;";
 414     out << "unset multiplot;";
 415
 416 }
 417
 418 void Calculator::filter(QMap<double, double> & map, const Doub & eps)
 419 {
 420     qDebug() << "filter crap smaller than " << eps;
 421     Int size, k;
 422     Doub diff, key;
 423     QList<double> values, list;
 424     list = map.uniqueKeys();
 425     //qDebug() << list;
 426     size = list.size();
 427     //qDebug() << "size of map:" << size;
 428     for(Int i = 0; i < size; i++) {
 429         key = list.at(i);
 430         //qDebug() << "["<<key<<"]" << "writing values into list";
 431         values = map.values(key);
 432         //qDebug() << "["<<key<<"]" << values.size() << "values for" << key;
 433         k = 0;
 434         for(Int k=0;k<values.size()-1;k++) {
 435             diff = abs(values.at(k+1) - values.at(k));
 436             //qDebug() << values.at(k+1) <<"-"<< values.at(k) << "=" << diff;
 437             if (diff <= eps) {
 438                 //qDebug() << diff << "<=" << eps;
 439                 //qDebug() << "["<<k<<"]" << "removing" << values.at(k) << "from list";
 440                 values.replace(k+1,values.at(k));
 441             }
 442         } // end for
 443         for (int j = 0; j < values.size(); j++)
 444             //qDebug() << "["<<j<<"]" << values.at(j);
 445         //qDebug() << "removing" << key << "from map";
 446         map.remove(key);
 447         for (int j = 0; j < values.size(); j++) {
 448             //qDebug() << "inserting" << values.at(j) << "into map at" << key;
 449             map.insertMulti(key,values.at(j));
 450         }
 451         //qDebug() << map.values(key);
 452         //qDebug() << endl << "end while" << endl;
 453     }
 454 }
 455
 456 void  Calculator::IC_setup(VecDoub &y, Doub x0, Doub y0, Doub xdot0, Doub ydot0,
 457                           Doub phi0, Doub params[])
 458 {
 459     //Int n = y.size();
 460     const Int sys = c->valuesMap.value(1);
 461
 462     /*
 463     qDebug() << "using " << "\tx0 = " << x0 << "\ty0 = " << y0
 464             << "\txdot0 = " << xdot0 << "\tydot0 = " << ydot0
 465             << "\tphi0 = " << fmod(phi0,360) << endl;
 466     */
 467
 468     switch(sys) {
 469         case 1:
 470         params[0] = c->valuesMap.value(2); // a;
 471         params[1] = c->valuesMap.value(3); // A;
 472         params[2] = c->valuesMap.value(4); // F;
 473         params[3] = c->valuesMap.value(5); // l;
 474         params[4] = c->valuesMap.value(6); // w;
 475         params[5] = c->valuesMap.value(25); // eta1;
 476         params[6] = 0;
 477         params[9] = c->valuesMap.value(32); // T;
 478         params[10] = 0;
 479         y[0] = x0, y[1] = xdot0;
 480         break;
 481         case 4:
 482         params[0] = c->valuesMap.value(2); // a;
 483         params[1] = fmod(c->valuesMap.value(3),360)*PI/180; // alpha;
 484         params[2] = c->valuesMap.value(4); // F;
 485         params[3] = c->valuesMap.value(5); // l;
 486         params[4] = c->valuesMap.value(6); // w;
 487         params[5] = c->valuesMap.value(25); // eta1;
 488         params[6] = c->valuesMap.value(26); // eta2;
 489         params[9] = c->valuesMap.value(32); // T;
 490         params[10] = fmod(c->valuesMap.value(40),360)*PI/180; // theta;
 491         params[17] = fmod(c->valuesMap.value(41),360)*PI/180; // psi;
 492         params[18] = c->valuesMap.value(42); // U;
 493         Doub e0 = 1.0/(params[5]+params[6]);
 494         Doub e1 = 1.0/(params[5]*params[3]);
 495         Doub e2 = 1.0/(params[6]*params[3]);
 496         Doub e3 = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 497         Doub fx = cos(params[10])*params[2];
 498         Doub fy = sin(params[10])*params[2];
 499         params[11] = e0; // 1/eta
 500         params[12] = e1;
 501         params[13] = e2;
 502         params[14] = e3;
 503         params[15] = fx;
 504         params[16] = fy;
 505         params[7] = params[3]*params[11]*params[5]; // l1;
 506         params[8] = params[3]*params[11]*params[6]; // l2;
 507         y[0] = x0, y[1] = y0, y[2] = phi0;
 508         break;
 509     }
 510     /*
 511     for (Int i = 0; i<n;i++)
 512         qDebug() << "y("<<i<<") =" << y[i];
 513     qDebug() << "eta1 = " << eta1 << "\teta2 = " << eta2 << "\tA = " << A << "\tF = " << F;
 514     */
 515 }
 516
 517 void Calculator::bifparam(const int &v, const QString &param, const Doub &r,
 518                           Doub params[])
 519 {
 520     switch(v) {
 521         case 1:
 522         if (param == "E")
 523             params[5] = r;
 524         else if (param == "F")
 525             params[2] = r;
 526         else if (param == "A")
 527             params[1] = r;
 528         else
 529             throw ("Fatal error!");
 530         break;
 531         case 4:
 532         if (param == "A")
 533             params[1] = fmod(r,360)*PI/180;
 534         else if (param == "F") {
 535             params[2] = r;
 536             params[15] = cos(params[10])*params[2];
 537             params[16] = sin(params[10])*params[2];
 538         }
 539         else if (param == "a")
 540             params[0] = r;
 541         else if (param == "w")
 542             params[4] = r;
 543         else if (param == "l") {
 544             params[3] = r;
 545             params[11] = 1.0/(params[5]+params[6]);
 546             params[12] = 1.0/(params[5]*params[3]);
 547             params[13] = 1.0/(params[6]*params[3]);
 548             params[14] = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 549             params[7] = params[3]*params[11]*params[5]; // l1;
 550             params[8] = params[3]*params[11]*params[6]; // l2;
 551         }
 552         else if (param == "e") {
 553             params[5] = r;
 554             params[11] = 1.0/(params[5]+params[6]);
 555             params[12] = 1.0/(params[5]*params[3]);
 556             params[13] = 1.0/(params[6]*params[3]);
 557             params[14] = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 558             params[7] = params[3]*params[11]*params[5]; // l1;
 559             params[8] = params[3]*params[11]*params[6]; // l2;
 560         }
 561         else if (param == "t") {
 562             params[10] = fmod(r,360)*PI/180;
 563             params[15] = cos(params[10])*params[2];
 564             params[16] = sin(params[10])*params[2];
 565         }
 566         else
 567             throw("Fatal error!");
 568         break;
 569         qDebug() << "Bifurcation parameter is " << param;
 570     }
 571 }
 572
 573 void Calculator::set_color(const Doub vel, QString & color)
 574 {
 575     Doub ctol = 0.05;
 576     if (vel > 1.0-ctol && vel < 1.0+ctol)
 577         color = "red";
 578     else if (vel > 0.5-ctol && vel < 0.5+ctol)
 579         color = "orange";
 580     else if (vel > 0.0-ctol && vel < 0.0+ctol)
 581         color = "grey";
 582     else if (vel > -0.5-ctol && vel < -0.5+ctol)
 583         color = "green";
 584     else if (vel > -1.0-ctol && vel < -1.0+ctol)
 585         color = "blue";
 586     else color = "black";
 587 }
 588
 589 void Calculator::calc2(const Int n)
 590 {
 591     // setting up a few start values
 592     qDebug() << "calc2 has been called: " << "calc2("<<n<<")";
 593     Normaldev gauss(0,1,time(NULL));
 594     Ran ran(time(NULL));
 595     Doub x1,x2,y1,y2,xtest;
 596     Doub L, tau, sigma;
 597     Doub x0, xdot0, y0, ydot0, phi0;
 598     x0 = c->valuesMap.value(8);
 599     xdot0 = c->valuesMap.value(10);
 600     ydot0 = c->valuesMap.value(24);
 601     y0 = c->valuesMap.value(23);
 602     phi0 = fmod(c->valuesMap.value(27),360)*PI/180;
 603     const Int samples = c->valuesMap.value(35);
 604     const Int v = c->valuesMap.value(1);
 605     const Int M = c->valuesMap.value(11);
 606     const Doub fromt = c->valuesMap.value(22);
 607     const Int s = c->valuesMap.value(14);
 608     const Doub h = c->valuesMap.value(13);
 609     const Doub t0 = c->valuesMap.value(7);
 610     const Doub fromr1 = c->valuesMap.value(17);
 611     const Doub fromr2 = c->valuesMap.value(37);
 612     const Doub tor1 = c->valuesMap.value(18);
 613     const Doub tor2 = c->valuesMap.value(38);
 614     const Doub stepr1 = c->valuesMap.value(19);
 615     const Doub stepr2 = c->valuesMap.value(39);
 616     const Doub eps = c->valuesMap.value(36);
 617     const Doub from = c->valuesMap.value(15);
 618     const Doub to = c->valuesMap.value(16);
 619     const Doub hstep = c->valuesMap.value(43);
 620     Doub params[256];
 621     //const Doub L = 1;
 622     const bool ruk = c->ruk;
 623     const bool eul = c->eul;
 624     QString param1 = c->str1;
 625     QString param2 = c->str2;
 626     qDebug() << "System" << v;
 627
 628     Doub x, hnew;
 629     Doub trajectory_length, deflection, velocity_angle, bias_angle;
 630     Int k, hit;
 631
 632     Doub xpi[samples][s];
 633     Doub ypi[samples][s];
 634     Doub phipi[samples][s];
 635     VecDoub y(n,0.0);
 636     VecDoub dydx(n,0.0);
 637     VecDoub dW(n, 0.0);
 638     VecDoub yout(n,0.0);
 639     VecDoub avg(n,0.0);
 640
 641     if (ruk) {
 642         switch(M) {
 643             case 0:
 644             qDebug() << "standard mode";
 645             qDebug() << "stepsize:" << h;
 646             qDebug() << "number of steps:" << s;
 647             x = t0;
 648             IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 649             qDebug() << phi0 << y[2];
 650             k = 0;
 651             L = 2*PI;
 652             tau = 2*PI/params[4];
 653             while(k < s) {
 654                 xtest = x+h;
 655                 if (xtest >= k*tau) {
 656                     k++;
 657                 }
 658                 derivs(x,y,dydx,params); // sets first dydx[] for rk4()
 659                 rk4(y,dydx,x,h,yout,params);
 660                 for (Int j=0;j<n;j++)
 661                     y[j] = yout[j];
 662                 x += h;
 663                 if (k >= s*fromt) {
 664                     switch(v) {
 665                         case 1:
 666                         rkValues_x.insert(x,yout[0]);
 667                         break;
 668                         case 4:
 669                         //xpi = ABS(fmod(yout[0],(2*PI)));
 670                         //ypi = ABS(fmod(yout[1],(2*PI)));
 671                         rkValues_x.insert(x,yout[0]);
 672                         rkValues_y.insert(x,yout[1]);
 673                         //rkValues_x.insert(x,xpi);
 674                         //rkValues_y.insert(x,ypi);
 675                         rkValues_phi.insert(x,(180/PI)*ABS(fmod(yout[2],(2*PI))));
 676                         x1 = y[0] + params[7]*cos(yout[2]);
 677                         y1 = y[1] + params[7]*sin(yout[2]);
 678                         x2 = y[0] - params[8]*cos(yout[2]);
 679                         y2 = y[1] - params[8]*sin(yout[2]);
 680                         rkValues_p1x.insert(x,x1);
 681                         rkValues_p1y.insert(x,y1);
 682                         rkValues_p2x.insert(x,x2);
 683                         rkValues_p2y.insert(x,y2);
 684                         rkValues_vx.insert(x,dydx[0]);
 685                         rkValues_vy.insert(x,dydx[1]);
 686                         rkValues_vphi.insert(x,dydx[2]);
 687                         break;
 688                     }
 689                 }
 690             } // end while
 691             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 692             qDebug() << "sizeof(rkValues_y) =" << rkValues_y.size();
 693             qDebug() << "sizeof(rkValues_phi) =" << rkValues_phi.size();
 694             write_to_file(file_rk_out, rkValues_x, rkValues_y, rkValues_phi, params);
 695             write_to_file(file_rk_pt1, rkValues_p1x, rkValues_p1y, rkValues_phi, params);
 696             write_to_file(file_rk_pt2, rkValues_p2x, rkValues_p2y, rkValues_phi, params);
 697             write_to_file(file_rk_vel, rkValues_vx, rkValues_vy, rkValues_vphi, params);
 698             break;
 699             case 1:
 700             qDebug() << "converge mode";
 701             Int iter;
 702             Doub hh, ss;
 703             for (iter=1;from+(iter*hstep)<=to;iter++) {
 704                 for (Int z=0;z<samples;z++) {
 705                     qDebug() << "at step" << z << "of" << samples;
 706                     x = t0;
 707                     hh = from+(iter*hstep);
 708                     ss = (s/hh);
 709                     qDebug() << "stepsize is" << hh << "with" << ss << "total steps.";
 710                     //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 711                     IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 712                     k = 0;
 713                     L = 2*PI;
 714                     tau = 2*PI/params[4];
 715                     while(k < s) {
 716                         //qDebug() << k << "<" << s;
 717                         //qDebug() << k << tau << k*tau;
 718                         xtest = x+hh;
 719                         if (xtest >= k*tau) {
 720                             //qDebug() << xtest << ">=" << k << "*" << tau << "("<<k*tau<<")";
 721                             hnew = k*tau-x;
 722                             derivs(x,y,dydx,params);
 723                             rk4(y,dydx,x,hnew,yout,params);
 724                             k++;
 725                             x += hnew;
 726                             //qDebug() << x;
 727                         } else {
 728                             //qDebug() << xtest << "=" << x<<"+"<<hh;
 729                             derivs(x,y,dydx,params);
 730                             rk4(y,dydx,x,hh,yout,params); // xold
 731                             x += hh;
 732                         }
 733                         for (Int j=0;j<n;j++)
 734                             y[j] = yout[j];
 735                     } // end while
 736                     switch(v) {
 737                         case 1:
 738                         rkValues_x.insert(hh,yout[0]);
 739                         break;
 740                         case 4:
 741                         rkValues_x.insert(hh,yout[0]);
 742                         rkValues_y.insert(hh,yout[1]);
 743                         rkValues_phi.insert(hh,yout[2]);
 744                         break;
 745                     }
 746                     x = t0;
 747                     //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 748                     IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 749                     k = 0;
 750                     L = 2*PI;
 751                     tau = 2*PI/params[4];
 752                     while (k < s) {
 753                         xtest = x+hh;
 754                         if (xtest >= k*tau) {
 755                             hnew = k*tau-x;
 756                             derivs(x,y,dydx,params);
 757                             derivw(x,y,dW,params,gauss);
 758                             euler(y,dydx,dW,x,hnew,yout);
 759                             k++;
 760                             x += hnew;
 761                         } else {
 762                             derivs(x,y,dydx,params);
 763                             derivw(x,y,dW,params,gauss);
 764                             euler(y,dydx,dW,x,hh,yout);
 765                             x += hh;
 766                         }
 767                         for (Int j=0;j<n;j++)
 768                             y[j] = yout[j];
 769                     } // end while
 770                     switch(v) {
 771                         case 1:
 772                         euValues_x.insert(hh,yout[0]);
 773                         break;
 774                         case 4:
 775                         euValues_x.insert(hh,yout[0]);
 776                         euValues_y.insert(hh,yout[1]);
 777                         euValues_phi.insert(hh,yout[2]);
 778                         break;
 779                     } // end switch
 780                 } // end samples
 781             } // end for
 782             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 783             write_to_file(file_rk_conv, rkValues_x, rkValues_y, rkValues_phi, params);
 784             write_to_file(file_eu_conv, euValues_x, euValues_y, euValues_phi, params);
 785             break;
 786             case 2:
 787             qDebug() << "bifurcation mode";
 788             qDebug() << "bifurcation parameter: " << param1;
 789             qDebug() << "Number of samples : " << samples;
 790             Doub r;
 791             for (Int j=0;(j*stepr1)<=(tor1-fromr1);j++) {
 792                 r = fromr1+j*stepr1;
 793                 qDebug() << "at step" << j << "of" << (tor1-fromr1)/stepr1;
 794                 for (Int z=0;z<samples;z++) {
 795                     x = t0;
 796                     x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 797                     //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
 798                     IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 799                     bifparam(v,param1,r,params);
 800                     //qDebug() << x << y[0] << y[1] << y[2];
 801                     //for(Int i = 0;i<=10;i++)
 802                       //  qDebug() << params[i];
 803                     // inner for: iterate over runge-kutta
 804                     k = 0;
 805                     L = 2*PI;
 806                     tau = 2*PI/params[4];
 807                     Doub hnew;
 808                     //qDebug() << "tau: " << tau << "\tL: " << L;
 809                     while (k<s) {
 810                         xtest = x+h; // test next x increment
 811                         //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
 812                         if (xtest >= k*tau) {
 813                             //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
 814                             //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
 815                             hnew = k*tau-x;
 816                             //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 817                             derivs(x,y,dydx,params);
 818                             rk4(y,dydx,x,hnew,yout,params);
 819                             switch(v) {
 820                                 case 1:
 821                                 xpi[z][k] = ABS(fmod(yout[0],(2*PI)));
 822                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 823                                 if (k >= s*fromt) {
 824                                     rkValues_x.insert(r,xpi[z][k]);
 825                                     rkValues_xdot.insert(r,dydx[0]);
 826                                     }
 827                                 break;
 828                                 case 4:
 829                                 xpi[z][k] = yout[0];
 830                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 831                                 ypi[z][k] = yout[1];
 832                                 phipi[z][k] = yout[2];
 833                                 rkDist_x.insert(r,SQR(yout[0]-x0));
 834                                 rkDist_y.insert(r,SQR(yout[1]-y0));
 835                                 rkDist_phi.insert(r,SQR(yout[2]-phi0));
 836                                 if (k >= s*fromt) {
 837                                     rkValues_x.insert(r,ABS(fmod(xpi[z][k],L)));
 838                                     rkValues_y.insert(r,ABS(fmod(ypi[z][k],L)));
 839                                     rkValues_phi.insert(r,ABS(fmod(phipi[z][k],L)));
 840                                     rkValues_xdot.insert(r,dydx[0]);
 841                                     rkValues_ydot.insert(r,dydx[1]);
 842                                     rkValues_phidot.insert(r,dydx[2]);
 843                                 }
 844                                 break;
 845                             } // end switch
 846                             k++;
 847                         }
 848                         derivs(x,y,dydx,params);
 849                         rk4(y,dydx,x,h,yout,params); // xold
 850                         for (Int j=0;j<n;j++)
 851                             y[j] = yout[j]; // y[x+h] !!!!!!!
 852                         x += h; // do it for real
 853                     } // end while
 854                     //qDebug() << "at" << x << "in time for k = " << k;
 855                     switch(v) {
 856                         Int m_x, m_y, m_phi;
 857                         case 1:
 858                         m_x = yout[0]/L;
 859                         avg[0] = (Doub) m_x/k;
 860                         rkAvgvelo_x.insert(r,avg[0]);
 861                         break;
 862                         case 4:
 863                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
 864                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
 865                         // makes only sense with r = params[10]
 866                         trajectory_length = sqrt(SQR(yout[0]-x0)+SQR(yout[1])-y0);
 867                         if (abs(yout[0]-x0) > trajectory_length) {
 868                             velocity_angle = PI;
 869                         } else {
 870                             velocity_angle = (acos((yout[0]-x0)/trajectory_length));
 871                         }
 872                         bias_angle = params[10];
 873                         deflection = (180/PI)*ABS(velocity_angle - bias_angle);
 874                         /*
 875                         qDebug() << velocity_angle << "("<<yout[0]-x0<<")"
 876                                 << "-" << bias_angle << "="
 877                                 << deflection
 878                                 << "[ length:"<<trajectory_length<<"]";
 879                         */
 880                         rkAvgvelo_x.insert(r,avg[0]);
 881                         rkAvgvelo_y.insert(r,avg[1]);
 882                         rkAvgvelo_phi.insert(r,avg[2]);
 883                         rkDefAngle.insert(r,deflection);
 884                         break;
 885                     } // end switch
 886                 } // samples loop
 887                 VecDoub msd_x(samples,0.0);
 888                 VecDoub msd_y(samples,0.0);
 889                 VecDoub msd_phi(samples,0.0);
 890                 Doub sigma_x[s-1];
 891                 Doub sigma_y[s-1];
 892                 Doub sigma_phi[s-1];
 893                 //qDebug() << "samples:" << samples << "steps:" << s;
 894                 for(Int i=0;i<s;i++) {
 895                     //qDebug() << "["<<i<<"]:---------------------------------";
 896                     for(Int j=0;j<samples;j++) {
 897                         //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
 898                         msd_x[j] = xpi[j][i];
 899                         //qDebug() << "["<<j<<"]:" << sper_x[j];
 900                         msd_y[j] = ypi[j][i];
 901                         msd_phi[j] = phipi[j][i];
 902                     }
 903                     diffusion(msd_x, samples, sigma);
 904                     sigma_x[i] = sigma;
 905                     diffusion(msd_y, samples, sigma);
 906                     sigma_y[i] = sigma;
 907                     diffusion(msd_phi, samples, sigma);
 908                     sigma_phi[i] = sigma;
 909                     //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
 910                 }
 911                 // the lists of keys in QMultiMap are read out backwards...
 912                 for(Int k=(s-1);k>=0;k--) {
 913                     rkVariance_x.insert(r,sigma_x[k]);
 914                     rkVariance_y.insert(r,sigma_y[k]);
 915                     rkVariance_phi.insert(r,sigma_phi[k]);
 916                 }
 917                //qDebug() << rkVariance_x;
 918             } // end outer for
 919             /*
 920              * we need to do something about this
 921              */
 922             /*
 923             filter(rkValues_x,eps);
 924             filter(rkValues_y,eps);
 925             filter(rkValues_xdot,eps);
 926             filter(rkValues_ydot,eps);
 927             filter(rkAvgvelo_x, eps);
 928             filter(rkAvgvelo_y, eps);
 929             filter(rkAvgvelo_phi,eps);*/
 930             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 931             qDebug() << "sizeof(rkValues_y) =" << rkValues_y.size();
 932             qDebug() << "sizeof(rkValues_xdot) =" << rkValues_xdot.size();
 933             qDebug() << "sizeof(rkValues_ydot) =" << rkValues_ydot.size();
 934             write_to_file(file_rk_bif, rkValues_x, rkValues_y, rkValues_phi, params);
 935             write_to_file(file_rk_biv, rkValues_xdot, rkValues_ydot, rkValues_phidot, params);
 936             write_to_file(file_rk_avg, rkAvgvelo_x, rkAvgvelo_y, rkAvgvelo_phi, params);
 937             write_to_file(file_rk_variance, rkVariance_x, rkVariance_y, rkVariance_phi,params);
 938             write_to_file(file_rk_distance, rkDist_x, rkDist_y, rkDist_phi, params);
 939             write_to_file(file_rk_defangle, rkDefAngle, NullMap, NullMap, params);
 940             break;
 941             case 3:
 942             qDebug() << "dual bifurcation mode";
 943             qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
 944             qDebug() << "Number of samples : " << samples;
 945             Doub stp1, stp2;
 946             QString vx_color, vy_color;
 947             QString name = "tmp.dat";
 948             QFile file(name);
 949             QTextStream out(&file);
 950             file.open(QIODevice::WriteOnly);
 951             out << "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl;
 952             for(Int it1=0;(it1*stepr1)<=(tor1-fromr1);it1++) {
 953                 for(Int it2=0;(it2*stepr2)<=(tor2-fromr2);it2++) {
 954                     stp1 = fromr1+it1*stepr1;
 955                     stp2 = fromr2+it2*stepr2;
 956                     qDebug() << stp1 << "("<<tor1<<")" << stp2 << "("<<tor2<<")";
 957                     for (Int z=0;z<samples;z++) {
 958                         x = t0;
 959                         x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 960                         IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 961                         bifparam(v,param1,stp1,params);
 962                         bifparam(v,param2,stp2,params);
 963                         /*
 964                         out << "# ";
 965                         for (Int i=0;i<10;i++)
 966                             out << "params["<<i<<"] = " <<  params[i] << " ";
 967                         out << endl;
 968                         */
 969                         // inner for: iterate over runge-kutta
 970                         k = 0;
 971                         L = 2*PI;
 972                         tau = 2*PI/params[4];
 973                         Doub hnew;
 974                         //qDebug() << "tau: " << tau << "\tL: " << L;
 975                         while (k<s) {
 976                             xtest = x+h; // test next x increment
 977                             //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
 978                             if (xtest >= k*tau) {
 979                                 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
 980                                 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
 981                                 hnew = k*tau-x;
 982                                 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 983                                 derivs(x,y,dydx,params);
 984                                 rk4(y,dydx,x,hnew,yout,params);
 985                                 k++;
 986                             }
 987                             derivs(x,y,dydx,params);
 988                             rk4(y,dydx,x,h,yout,params); // xold
 989                             for (Int j=0;j<n;j++)
 990                                 y[j] = yout[j]; // y[x+h] !!!!!!!
 991                             x += h; // do it for real
 992                         } // end while
 993                         Int m_x, m_y, m_phi;
 994                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
 995                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k;
 996                         avg[2] = (Doub) m_phi/k;
 997                         set_color(avg[0], vx_color);
 998                         set_color(avg[1], vy_color);
 999                         out << stp1 << "\t" << stp2 << "\t"
1000                                 << avg[0] <<  "\t" << vx_color << "\t" << avg[1]
1001                                 << "\t" << vy_color << "\t" << avg[2] << endl;
1002                     } // end samples
1003                 } // end 2. outer for
1004             } // end 1. outer for
1005             file.close();
1006             qDebug() << "moving" << name << "to safety";
1007             file.copy(name, file_rk_dualbif);
1008             break;
1009         } // end switch
1010     } // end if
1011
1012     if (eul) {
1013         switch(M) {
1014           case 0:
1015           qDebug() << "(euler) standard mode";
1016           qDebug() << "stepsize:" << h;
1017           qDebug() << "number of steps:" << s;
1018           x = t0;
1019           IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1020           k = 0;
1021           hit = 0;
1022           L = 2*PI;
1023           tau = 2*PI/params[4];
1024           while(k < s) {
1025               //qDebug() << k << "<" << s;
1026               //qDebug() << "x is now" << x;
1027               //qDebug() << "deriving values at" << x;
1028               //qDebug() << x0, y0, phi0;
1029               derivs(x,y,dydx,params); // sets first dydx[] for rk4()
1030               derivw(x,y,dW,params,gauss);
1031               euler(y,dydx,dW,x,h,yout);
1032               //for (Int i=0;i<=18;i++)
1033                 //qDebug() << "params["<<i<<"]" << params[i];
1034               if (x >= hit*0.1) {
1035                   //qDebug() << x <<"=="<< hit*0.1;
1036                   //qDebug() << "writing2 to map";
1037                   hit++;
1038                   //qDebug() << "next value to write out" << hit*0.1;
1039                   if (k >= s*fromt) {
1040                       switch(v) {
1041                       case 1:
1042                           euValues_x.insert(x,yout[0]);
1043                       break;
1044                       case 4:
1045                         //xpi = ABS(fmod(yout[0],(2*PI)));
1046                         //ypi = ABS(fmod(yout[1],(2*PI)));
1047                         euValues_x.insert(x,yout[0]);
1048                         euValues_y.insert(x,yout[1]);
1049                         //rkValues_x.insert(x,xpi);
1050                         //rkValues_y.insert(x,ypi);
1051                         euValues_phi.insert(x,yout[2]);
1052                         x1 = yout[0] + params[7]*cos(yout[2]);
1053                         y1 = yout[1] + params[7]*sin(yout[2]);
1054                         x2 = yout[0] - params[8]*cos(yout[2]);
1055                         y2 = yout[1] - params[8]*sin(yout[2]);
1056                         euValues_p1x.insert(x,x1);
1057                         euValues_p1y.insert(x,y1);
1058                         euValues_p2x.insert(x,x2);
1059                         euValues_p2y.insert(x,y2);
1060                         euValues_vx.insert(x,dydx[0]);
1061                         euValues_vy.insert(x,dydx[1]);
1062                         euValues_vphi.insert(x,dydx[2]);
1063                         break;
1064                     } // end fromt switch
1065                   } // end fromt if
1066               } // end hit if
1067               //qDebug() << "making step of size" << h;
1068               for (Int j=0;j<n;j++)
1069                   y[j] = yout[j];
1070               x += h;
1071               //qDebug() << "new x is now" << x;
1072               if (x >= k*tau) {
1073                   //qDebug() << x << ">=" << k*tau;
1074                   //qDebug() << "incrementing k";
1075                   k++;
1076                   //qDebug() << "new value of x needed is" << k*tau;
1077               }
1078           } // end while
1079           qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
1080           qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
1081           qDebug() << "sizeof(euValues_phi) =" << euValues_phi.size();
1082           write_to_file(file_eu_out, euValues_x, euValues_y, euValues_phi, params);
1083           write_to_file(file_eu_pt1, euValues_p1x, euValues_p1y, euValues_phi, params);
1084           write_to_file(file_eu_pt2, euValues_p2x, euValues_p2y, euValues_phi, params);
1085           write_to_file(file_eu_vel, euValues_vx, euValues_vy, euValues_vphi, params);
1086           break;
1087           case 2:
1088           qDebug() << "bifurcation mode";
1089           qDebug() << "bifurcation parameter: " << param1;
1090           qDebug() << "Number of samples : " << samples;
1091           Doub r;
1092           for (Int j=0;(j*stepr1)<=(tor1-fromr1);j++) {
1093               r = fromr1+j*stepr1;
1094               qDebug() << "at step" << j << "of" << (tor1-fromr1)/stepr1;
1095               for (Int z=0;z<samples;z++) {
1096                   x = t0;
1097                   x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
1098                   //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
1099                   IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1100                   bifparam(v,param1,r,params);
1101                   // inner for: iterate over runge-kutta
1102                   k = 0;
1103                   L = 2*PI;
1104                   tau = 2*PI/params[4];
1105                   Doub hnew;
1106                   //qDebug() << "tau: " << tau << "\tL: " << L;
1107                   while (k<s) {
1108                       xtest = x+h; // test next x increment
1109                       //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1110                       if (xtest >= k*tau) {
1111                           //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1112                           //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1113                           hnew = k*tau-x;
1114                           //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1115                           derivs(x,y,dydx,params);
1116                           derivw(x,y, dW, params,gauss);
1117                           euler(y,dydx,dW,x,hnew,yout);
1118                           switch(v) {
1119                           case 1:
1120                               xpi[z][k] = yout[0];
1121                               if (k >= s*fromt) {
1122                                   euValues_x.insert(r,ABS(fmod(xpi[z][k],L)));
1123                                   euValues_xdot.insert(r,dydx[0]);
1124                               }
1125                               break;
1126                           case 4:
1127                               xpi[z][k] = yout[0];
1128                               //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
1129                               ypi[z][k] = yout[1];
1130                               phipi[z][k] = yout[2];
1131                               euDist_x.insert(r,SQR(yout[0]-x0));
1132                               euDist_y.insert(r,SQR(yout[1]-y0));
1133                               euDist_phi.insert(r,SQR(yout[2]-phi0));
1134                               if (k >= s*fromt) {
1135                                   euValues_x.insert(r,ABS(fmod(xpi[z][k],L)));
1136                                   euValues_y.insert(r,ABS(fmod(ypi[z][k],L)));
1137                                   euValues_phi.insert(r,ABS(fmod(phipi[z][k],L)));
1138                                   euValues_xdot.insert(r,dydx[0]);
1139                                   euValues_ydot.insert(r,dydx[1]);
1140                                   euValues_phidot.insert(r,dydx[2]);
1141                               }
1142                               break;
1143                           } // end switch
1144                           k++;
1145                       }
1146                       derivs(x,y,dydx,params);
1147                       derivw(x,y,dW, params,gauss);
1148                       euler(y,dydx,dW,x,h,yout); // xold
1149                       for (Int j=0;j<n;j++)
1150                           y[j] = yout[j]; // y[x+h] !!!!!!!
1151                       x += h; // do it for real
1152                   } // end while
1153                   //qDebug() << "at" << x << "in time for k = " << k;
1154                   switch(v) {
1155                       Int m_x, m_y, m_phi;
1156                       case 1:
1157                       m_x = yout[0]/L;
1158                       avg[0] = (Doub) m_x/k;
1159                       euAvgvelo_x.insert(r,avg[0]);
1160                       break;
1161                       case 4:
1162                       m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
1163                       avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
1164                       trajectory_length = sqrt(SQR(yout[0]-x0)+SQR(yout[1])-y0);
1165                       if (abs(yout[0]-x0) > trajectory_length) {
1166                           velocity_angle = PI;
1167                       } else {
1168                           velocity_angle = (acos((yout[0]-x0)/trajectory_length));
1169                       }
1170                       bias_angle = params[10];
1171                       deflection = (180/PI)*ABS(velocity_angle - bias_angle);
1172                       euAvgvelo_x.insert(r,avg[0]);
1173                       euAvgvelo_y.insert(r,avg[1]);
1174                       euAvgvelo_phi.insert(r,avg[2]);
1175                       euDefAngle.insert(r,deflection);
1176                       break;
1177                   } // end switch
1178               } // samples loop
1179               VecDoub msd_x(samples,0.0);
1180               VecDoub msd_y(samples,0.0);
1181               VecDoub msd_phi(samples,0.0);
1182               Doub sigma_x[s-1];
1183               Doub sigma_y[s-1];
1184               Doub sigma_phi[s-1];
1185               //qDebug() << "samples:" << samples << "steps:" << s;
1186               for(Int i=0;i<s;i++) {
1187                   //qDebug() << "["<<i<<"]:---------------------------------";
1188                   for(Int j=0;j<samples;j++) {
1189                       //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
1190                       msd_x[j] = xpi[j][i];
1191                       //qDebug() << "["<<j<<"]:" << sper_x[j];
1192                       msd_y[j] = ypi[j][i];
1193                       msd_phi[j] = phipi[j][i];
1194                   }
1195                   diffusion(msd_x, samples, sigma);
1196                   sigma_x[i] = sigma;
1197                   diffusion(msd_y, samples, sigma);
1198                   sigma_y[i] = sigma;
1199                   diffusion(msd_phi, samples, sigma);
1200                   sigma_phi[i] = sigma;
1201                   //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
1202               }
1203               // the lists of keys in QMultiMap are read out backwards...
1204               for(Int k=(s-1);k>=0;k--) {
1205                   euVariance_x.insert(r,sigma_x[k]);
1206                   euVariance_y.insert(r,sigma_y[k]);
1207                   euVariance_phi.insert(r,sigma_phi[k]);
1208               }
1209               //qDebug() << euVariance_x;
1210           } // end outer for
1211           /*
1212            * we need to do something about this
1213            */
1214           qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
1215           qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
1216           qDebug() << "sizeof(euValues_xdot) =" << euValues_xdot.size();
1217           qDebug() << "sizeof(euValues_ydot) =" << euValues_ydot.size();
1218           write_to_file(file_eu_bif, euValues_x, euValues_y, euValues_phi, params);
1219           write_to_file(file_eu_biv, euValues_xdot, euValues_ydot, euValues_phidot, params);
1220           write_to_file(file_eu_avg, euAvgvelo_x, euAvgvelo_y, euAvgvelo_phi, params);
1221           write_to_file(file_eu_variance, euVariance_x, euVariance_y, euVariance_phi,params);
1222           write_to_file(file_eu_distance, euDist_x, euDist_y, euDist_phi, params);
1223           write_to_file(file_eu_defangle, euDefAngle, NullMap, NullMap, params);
1224           break;
1225           case 3:
1226           qDebug() << "dual bifurcation mode";
1227           qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
1228           qDebug() << "Number of samples : " << samples;
1229           Doub stp1,stp2;
1230           QString vx_color, vy_color;
1231           QString name = "tmp.dat";
1232           QFile file(name);
1233           QTextStream out(&file);
1234           file.open(QIODevice::WriteOnly);
1235           out << "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl;
1236           for (Int it1=0;(it1*stepr1)<=(tor1-fromr1);it1++) {
1237               for (Int it2=0;(it2*stepr2)<=(tor2-fromr2);it2++) {
1238                   stp1 = fromr1+it1*stepr1;
1239                   stp2 = fromr2+it2*stepr2;
1240                   qDebug() << stp1 << "("<<tor1<<")" << stp2 << "("<<tor2<<")";
1241                   for (Int z=0;z<samples;z++) {
1242                       x = t0;
1243                       x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
1244                       IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1245                       bifparam(v,param1,stp1,params);
1246                       bifparam(v,param2,stp2,params);
1247                       // inner for: iterate over runge-kutta
1248                       Int k = 0;
1249                       L = 2*PI;
1250                       tau = 2*PI/params[4];
1251                       Doub hnew;
1252                       //qDebug() << "tau: " << tau << "\tL: " << L;
1253                       while (k<s) {
1254                           xtest = x+h; // test next x increment
1255                           //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1256                           if (xtest >= k*tau) {
1257                               //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1258                               //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1259                               hnew = k*tau-x;
1260                               //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1261                               derivs(x,y,dydx,params);
1262                               derivw(x,y, dW,params,gauss);
1263                               euler(y,dydx,dW,x,hnew,yout);
1264                               k++;
1265                           }
1266                           derivs(x,y,dydx,params);
1267                           derivw(x,y, dW,params,gauss);
1268                           euler(y,dydx,dW,x,h,yout); // xold
1269                           for (Int j=0;j<n;j++)
1270                               y[j] = yout[j]; // y[x+h] !!!!!!!
1271                           x += h; // do it for real
1272                       } // end while
1273                         Int m_x, m_y, m_phi;
1274                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
1275                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k;
1276                         avg[2] = (Doub) m_phi/k;
1277                         set_color(avg[0], vx_color);
1278                         set_color(avg[1], vy_color);
1279                         out << stp1 << "\t" << stp2 << "\t"
1280                                 << avg[0] <<  "\t" << vx_color << "\t" << avg[1]
1281                                 << "\t" << vy_color << "\t" << avg[2] << endl;
1282
1283                   } // end samples
1284               } // end 2. outer for
1285           } // end 1. outer for
1286           file.close();
1287           qDebug() << "moving" << name << "to safety";
1288           file.copy(name, file_rk_dualbif);
1289           break;
1290         } // end switch
1291     } // end if
1292
1293     if (!ruk && !eul)
1294         throw("No method specified!");
1295
1296     rkValues_x.clear();
1297     rkValues_xdot.clear();
1298     rkValues_y.clear();
1299     rkValues_ydot.clear();
1300     rkValues_phi.clear();
1301     rkValues_phidot.clear();
1302     rkAvgvelo_x.clear();
1303     rkAvgvelo_y.clear();
1304     rkAvgvelo_phi.clear();
1305
1306     euValues_x.clear();
1307     euValues_xdot.clear();
1308     euValues_y.clear();
1309     euValues_ydot.clear();
1310     euValues_phi.clear();
1311     euValues_phidot.clear();
1312     euAvgvelo_x.clear();
1313     euAvgvelo_y.clear();
1314     euAvgvelo_phi.clear();
1315
1316 } // end calculator
1317
1318 //=============================================================================