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 if (param == "T") {
 567             params[9] = r;
 568         }
 569         else
 570             throw("Fatal error!");
 571         break;
 572         qDebug() << "Bifurcation parameter is " << param;
 573     }
 574 }
 575
 576 void Calculator::set_color(const Doub vel, QString & color)
 577 {
 578     Doub ctol = 0.1;
 579     if (vel > 1.0-ctol && vel < 1.0+ctol)
 580         color = "red";
 581     else if (vel > 0.5-ctol && vel < 0.5+ctol)
 582         color = "orange";
 583     else if (vel > 0.0-ctol && vel < 0.0+ctol)
 584         color = "grey";
 585     else if (vel > -0.5-ctol && vel < -0.5+ctol)
 586         color = "green";
 587     else if (vel > -1.0-ctol && vel < -1.0+ctol)
 588         color = "blue";
 589     else color = "black";
 590 }
 591
 592 void Calculator::calc2(const Int n)
 593 {
 594     // setting up a few start values
 595     qDebug() << "calc2 has been called: " << "calc2("<<n<<")";
 596     Normaldev gauss(0,1,time(NULL));
 597     Ran ran(time(NULL));
 598     Doub x1,x2,y1,y2,xtest;
 599     Doub L, tau, sigma;
 600     Doub x0, xdot0, y0, ydot0, phi0;
 601     x0 = c->valuesMap.value(8);
 602     xdot0 = c->valuesMap.value(10);
 603     ydot0 = c->valuesMap.value(24);
 604     y0 = c->valuesMap.value(23);
 605     phi0 = fmod(c->valuesMap.value(27),360)*PI/180;
 606     const Int samples = c->valuesMap.value(35);
 607     const Int v = c->valuesMap.value(1);
 608     const Int M = c->valuesMap.value(11);
 609     const Doub fromt = c->valuesMap.value(22);
 610     const Int s = c->valuesMap.value(14);
 611     const Doub h = c->valuesMap.value(13);
 612     const Doub t0 = c->valuesMap.value(7);
 613     const Doub fromr1 = c->valuesMap.value(17);
 614     const Doub fromr2 = c->valuesMap.value(37);
 615     const Doub tor1 = c->valuesMap.value(18);
 616     const Doub tor2 = c->valuesMap.value(38);
 617     const Doub stepr1 = c->valuesMap.value(19);
 618     const Doub stepr2 = c->valuesMap.value(39);
 619     const Doub eps = c->valuesMap.value(36);
 620     const Doub from = c->valuesMap.value(15);
 621     const Doub to = c->valuesMap.value(16);
 622     const Doub hstep = c->valuesMap.value(43);
 623     Doub params[256];
 624     //const Doub L = 1;
 625     const bool ruk = c->ruk;
 626     const bool eul = c->eul;
 627     QString param1 = c->str1;
 628     QString param2 = c->str2;
 629     qDebug() << "System" << v;
 630
 631     Doub x, hnew;
 632     Doub trajectory_length, deflection, velocity_angle, bias_angle;
 633     Int k, hit;
 634
 635     Doub xpi[samples][s];
 636     Doub ypi[samples][s];
 637     Doub phipi[samples][s];
 638     VecDoub y(n,0.0);
 639     VecDoub dydx(n,0.0);
 640     VecDoub dW(n, 0.0);
 641     VecDoub yout(n,0.0);
 642     VecDoub avg(n,0.0);
 643
 644     if (ruk) {
 645         switch(M) {
 646             case 0:
 647             qDebug() << "standard mode";
 648             qDebug() << "stepsize:" << h;
 649             qDebug() << "number of steps:" << s;
 650             x = t0;
 651             IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 652             qDebug() << phi0 << y[2];
 653             k = 0;
 654             L = 2*PI;
 655             tau = 2*PI/params[4];
 656             while(k < s) {
 657                 xtest = x+h;
 658                 if (xtest >= k*tau) {
 659                     k++;
 660                 }
 661                 derivs(x,y,dydx,params); // sets first dydx[] for rk4()
 662                 rk4(y,dydx,x,h,yout,params);
 663                 for (Int j=0;j<n;j++)
 664                     y[j] = yout[j];
 665                 x += h;
 666                 if (k >= s*fromt) {
 667                     switch(v) {
 668                         case 1:
 669                         rkValues_x.insert(x,yout[0]);
 670                         break;
 671                         case 4:
 672                         //xpi = ABS(fmod(yout[0],(2*PI)));
 673                         //ypi = ABS(fmod(yout[1],(2*PI)));
 674                         rkValues_x.insert(x,yout[0]);
 675                         rkValues_y.insert(x,yout[1]);
 676                         //rkValues_x.insert(x,xpi);
 677                         //rkValues_y.insert(x,ypi);
 678                         rkValues_phi.insert(x,(180/PI)*ABS(fmod(yout[2],(2*PI))));
 679                         x1 = y[0] + params[7]*cos(yout[2]);
 680                         y1 = y[1] + params[7]*sin(yout[2]);
 681                         x2 = y[0] - params[8]*cos(yout[2]);
 682                         y2 = y[1] - params[8]*sin(yout[2]);
 683                         rkValues_p1x.insert(x,x1);
 684                         rkValues_p1y.insert(x,y1);
 685                         rkValues_p2x.insert(x,x2);
 686                         rkValues_p2y.insert(x,y2);
 687                         rkValues_vx.insert(x,dydx[0]);
 688                         rkValues_vy.insert(x,dydx[1]);
 689                         rkValues_vphi.insert(x,dydx[2]);
 690                         break;
 691                     }
 692                 }
 693             } // end while
 694             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 695             qDebug() << "sizeof(rkValues_y) =" << rkValues_y.size();
 696             qDebug() << "sizeof(rkValues_phi) =" << rkValues_phi.size();
 697             write_to_file(file_rk_out, rkValues_x, rkValues_y, rkValues_phi, params);
 698             write_to_file(file_rk_pt1, rkValues_p1x, rkValues_p1y, rkValues_phi, params);
 699             write_to_file(file_rk_pt2, rkValues_p2x, rkValues_p2y, rkValues_phi, params);
 700             write_to_file(file_rk_vel, rkValues_vx, rkValues_vy, rkValues_vphi, params);
 701             break;
 702             case 1:
 703             qDebug() << "converge mode";
 704             Int iter;
 705             Doub hh, ss;
 706             for (iter=1;from+(iter*hstep)<=to;iter++) {
 707                 for (Int z=0;z<samples;z++) {
 708                     qDebug() << "at step" << z << "of" << samples;
 709                     x = t0;
 710                     hh = from+(iter*hstep);
 711                     ss = (s/hh);
 712                     qDebug() << "stepsize is" << hh << "with" << ss << "total steps.";
 713                     //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 714                     IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 715                     k = 0;
 716                     L = 2*PI;
 717                     tau = 2*PI/params[4];
 718                     while(k < s) {
 719                         //qDebug() << k << "<" << s;
 720                         //qDebug() << k << tau << k*tau;
 721                         xtest = x+hh;
 722                         if (xtest >= k*tau) {
 723                             //qDebug() << xtest << ">=" << k << "*" << tau << "("<<k*tau<<")";
 724                             hnew = k*tau-x;
 725                             derivs(x,y,dydx,params);
 726                             rk4(y,dydx,x,hnew,yout,params);
 727                             k++;
 728                             x += hnew;
 729                             //qDebug() << x;
 730                         } else {
 731                             //qDebug() << xtest << "=" << x<<"+"<<hh;
 732                             derivs(x,y,dydx,params);
 733                             rk4(y,dydx,x,hh,yout,params); // xold
 734                             x += hh;
 735                         }
 736                         for (Int j=0;j<n;j++)
 737                             y[j] = yout[j];
 738                     } // end while
 739                     switch(v) {
 740                         case 1:
 741                         rkValues_x.insert(hh,yout[0]);
 742                         break;
 743                         case 4:
 744                         rkValues_x.insert(hh,yout[0]);
 745                         rkValues_y.insert(hh,yout[1]);
 746                         rkValues_phi.insert(hh,yout[2]);
 747                         break;
 748                     }
 749                     x = t0;
 750                     //x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 751                     IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 752                     k = 0;
 753                     L = 2*PI;
 754                     tau = 2*PI/params[4];
 755                     while (k < s) {
 756                         xtest = x+hh;
 757                         if (xtest >= k*tau) {
 758                             hnew = k*tau-x;
 759                             derivs(x,y,dydx,params);
 760                             derivw(x,y,dW,params,gauss);
 761                             euler(y,dydx,dW,x,hnew,yout);
 762                             k++;
 763                             x += hnew;
 764                         } else {
 765                             derivs(x,y,dydx,params);
 766                             derivw(x,y,dW,params,gauss);
 767                             euler(y,dydx,dW,x,hh,yout);
 768                             x += hh;
 769                         }
 770                         for (Int j=0;j<n;j++)
 771                             y[j] = yout[j];
 772                     } // end while
 773                     switch(v) {
 774                         case 1:
 775                         euValues_x.insert(hh,yout[0]);
 776                         break;
 777                         case 4:
 778                         euValues_x.insert(hh,yout[0]);
 779                         euValues_y.insert(hh,yout[1]);
 780                         euValues_phi.insert(hh,yout[2]);
 781                         break;
 782                     } // end switch
 783                 } // end samples
 784             } // end for
 785             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 786             write_to_file(file_rk_conv, rkValues_x, rkValues_y, rkValues_phi, params);
 787             write_to_file(file_eu_conv, euValues_x, euValues_y, euValues_phi, params);
 788             break;
 789             case 2:
 790             qDebug() << "bifurcation mode";
 791             qDebug() << "bifurcation parameter: " << param1;
 792             qDebug() << "Number of samples : " << samples;
 793             Doub r;
 794             for (Int j=0;(j*stepr1)<=(tor1-fromr1);j++) {
 795                 r = fromr1+j*stepr1;
 796                 qDebug() << "at step" << j << "of" << (tor1-fromr1)/stepr1;
 797                 for (Int z=0;z<samples;z++) {
 798                     x = t0;
 799                     x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 800                     //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
 801                     IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 802                     bifparam(v,param1,r,params);
 803                     //qDebug() << x << y[0] << y[1] << y[2];
 804                     //for(Int i = 0;i<=10;i++)
 805                       //  qDebug() << params[i];
 806                     // inner for: iterate over runge-kutta
 807                     k = 0;
 808                     L = 2*PI;
 809                     tau = 2*PI/params[4];
 810                     Doub hnew;
 811                     //qDebug() << "tau: " << tau << "\tL: " << L;
 812                     while (k<s) {
 813                         xtest = x+h; // test next x increment
 814                         //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
 815                         if (xtest >= k*tau) {
 816                             //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
 817                             //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
 818                             hnew = k*tau-x;
 819                             //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 820                             derivs(x,y,dydx,params);
 821                             rk4(y,dydx,x,hnew,yout,params);
 822                             switch(v) {
 823                                 case 1:
 824                                 xpi[z][k] = ABS(fmod(yout[0],(2*PI)));
 825                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 826                                 if (k >= s*fromt) {
 827                                     rkValues_x.insert(r,xpi[z][k]);
 828                                     rkValues_xdot.insert(r,dydx[0]);
 829                                     }
 830                                 break;
 831                                 case 4:
 832                                 xpi[z][k] = yout[0];
 833                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 834                                 ypi[z][k] = yout[1];
 835                                 phipi[z][k] = yout[2];
 836                                 rkDist_x.insert(r,SQR(yout[0]-x0));
 837                                 rkDist_y.insert(r,SQR(yout[1]-y0));
 838                                 rkDist_phi.insert(r,SQR(yout[2]-phi0));
 839                                 if (k >= s*fromt) {
 840                                     rkValues_x.insert(r,ABS(fmod(xpi[z][k],L)));
 841                                     rkValues_y.insert(r,ABS(fmod(ypi[z][k],L)));
 842                                     rkValues_phi.insert(r,ABS(fmod(phipi[z][k],L)));
 843                                     rkValues_xdot.insert(r,dydx[0]);
 844                                     rkValues_ydot.insert(r,dydx[1]);
 845                                     rkValues_phidot.insert(r,dydx[2]);
 846                                 }
 847                                 break;
 848                             } // end switch
 849                             k++;
 850                         }
 851                         derivs(x,y,dydx,params);
 852                         rk4(y,dydx,x,h,yout,params); // xold
 853                         for (Int j=0;j<n;j++)
 854                             y[j] = yout[j]; // y[x+h] !!!!!!!
 855                         x += h; // do it for real
 856                     } // end while
 857                     //qDebug() << "at" << x << "in time for k = " << k;
 858                     switch(v) {
 859                         Int m_x, m_y, m_phi;
 860                         case 1:
 861                         m_x = yout[0]/L;
 862                         avg[0] = (Doub) m_x/k;
 863                         rkAvgvelo_x.insert(r,avg[0]);
 864                         break;
 865                         case 4:
 866                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
 867                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
 868                         // makes only sense with r = params[10]
 869                         trajectory_length = sqrt(SQR(yout[0]-x0)+SQR(yout[1])-y0);
 870                         if (abs(yout[0]-x0) > trajectory_length) {
 871                             velocity_angle = PI;
 872                         } else {
 873                             velocity_angle = (acos((yout[0]-x0)/trajectory_length));
 874                         }
 875                         bias_angle = params[10];
 876                         deflection = (180/PI)*ABS(velocity_angle - bias_angle);
 877                         /*
 878                         qDebug() << velocity_angle << "("<<yout[0]-x0<<")"
 879                                 << "-" << bias_angle << "="
 880                                 << deflection
 881                                 << "[ length:"<<trajectory_length<<"]";
 882                         */
 883                         rkAvgvelo_x.insert(r,avg[0]);
 884                         rkAvgvelo_y.insert(r,avg[1]);
 885                         rkAvgvelo_phi.insert(r,avg[2]);
 886                         rkDefAngle.insert(r,deflection);
 887                         break;
 888                     } // end switch
 889                 } // samples loop
 890                 VecDoub msd_x(samples,0.0);
 891                 VecDoub msd_y(samples,0.0);
 892                 VecDoub msd_phi(samples,0.0);
 893                 Doub sigma_x[s-1];
 894                 Doub sigma_y[s-1];
 895                 Doub sigma_phi[s-1];
 896                 //qDebug() << "samples:" << samples << "steps:" << s;
 897                 for(Int i=0;i<s;i++) {
 898                     //qDebug() << "["<<i<<"]:---------------------------------";
 899                     for(Int j=0;j<samples;j++) {
 900                         //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
 901                         msd_x[j] = xpi[j][i];
 902                         //qDebug() << "["<<j<<"]:" << sper_x[j];
 903                         msd_y[j] = ypi[j][i];
 904                         msd_phi[j] = phipi[j][i];
 905                     }
 906                     diffusion(msd_x, samples, sigma);
 907                     sigma_x[i] = sigma;
 908                     diffusion(msd_y, samples, sigma);
 909                     sigma_y[i] = sigma;
 910                     diffusion(msd_phi, samples, sigma);
 911                     sigma_phi[i] = sigma;
 912                     //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
 913                 }
 914                 // the lists of keys in QMultiMap are read out backwards...
 915                 for(Int k=(s-1);k>=0;k--) {
 916                     rkVariance_x.insert(r,sigma_x[k]);
 917                     rkVariance_y.insert(r,sigma_y[k]);
 918                     rkVariance_phi.insert(r,sigma_phi[k]);
 919                 }
 920                //qDebug() << rkVariance_x;
 921             } // end outer for
 922             /*
 923              * we need to do something about this
 924              */
 925             /*
 926             filter(rkValues_x,eps);
 927             filter(rkValues_y,eps);
 928             filter(rkValues_xdot,eps);
 929             filter(rkValues_ydot,eps);
 930             filter(rkAvgvelo_x, eps);
 931             filter(rkAvgvelo_y, eps);
 932             filter(rkAvgvelo_phi,eps);*/
 933             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 934             qDebug() << "sizeof(rkValues_y) =" << rkValues_y.size();
 935             qDebug() << "sizeof(rkValues_xdot) =" << rkValues_xdot.size();
 936             qDebug() << "sizeof(rkValues_ydot) =" << rkValues_ydot.size();
 937             write_to_file(file_rk_bif, rkValues_x, rkValues_y, rkValues_phi, params);
 938             write_to_file(file_rk_biv, rkValues_xdot, rkValues_ydot, rkValues_phidot, params);
 939             write_to_file(file_rk_avg, rkAvgvelo_x, rkAvgvelo_y, rkAvgvelo_phi, params);
 940             write_to_file(file_rk_variance, rkVariance_x, rkVariance_y, rkVariance_phi,params);
 941             write_to_file(file_rk_distance, rkDist_x, rkDist_y, rkDist_phi, params);
 942             write_to_file(file_rk_defangle, rkDefAngle, NullMap, NullMap, params);
 943             break;
 944             case 3:
 945             qDebug() << "dual bifurcation mode";
 946             qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
 947             qDebug() << "Number of samples : " << samples;
 948             Doub stp1, stp2;
 949             QString vx_color, vy_color;
 950             QString name = "tmp.dat";
 951             QFile file(name);
 952             QTextStream out(&file);
 953             file.open(QIODevice::WriteOnly);
 954             out << "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl;
 955             for(Int it1=0;(it1*stepr1)<=(tor1-fromr1);it1++) {
 956                 for(Int it2=0;(it2*stepr2)<=(tor2-fromr2);it2++) {
 957                     stp1 = fromr1+it1*stepr1;
 958                     stp2 = fromr2+it2*stepr2;
 959                     qDebug() << stp1 << "("<<tor1<<")" << stp2 << "("<<tor2<<")";
 960                     for (Int z=0;z<samples;z++) {
 961                         x = t0;
 962                         x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 963                         IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 964                         bifparam(v,param1,stp1,params);
 965                         bifparam(v,param2,stp2,params);
 966                         /*
 967                         out << "# ";
 968                         for (Int i=0;i<10;i++)
 969                             out << "params["<<i<<"] = " <<  params[i] << " ";
 970                         out << endl;
 971                         */
 972                         // inner for: iterate over runge-kutta
 973                         k = 0;
 974                         L = 2*PI;
 975                         tau = 2*PI/params[4];
 976                         Doub hnew;
 977                         //qDebug() << "tau: " << tau << "\tL: " << L;
 978                         while (k<s) {
 979                             xtest = x+h; // test next x increment
 980                             //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
 981                             if (xtest >= k*tau) {
 982                                 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
 983                                 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
 984                                 hnew = k*tau-x;
 985                                 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 986                                 derivs(x,y,dydx,params);
 987                                 rk4(y,dydx,x,hnew,yout,params);
 988                                 k++;
 989                             }
 990                             derivs(x,y,dydx,params);
 991                             rk4(y,dydx,x,h,yout,params); // xold
 992                             for (Int j=0;j<n;j++)
 993                                 y[j] = yout[j]; // y[x+h] !!!!!!!
 994                             x += h; // do it for real
 995                         } // end while
 996                         Int m_x, m_y, m_phi;
 997                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
 998                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k;
 999                         avg[2] = (Doub) m_phi/k;
1000                         set_color(avg[0], vx_color);
1001                         set_color(avg[1], vy_color);
1002                         out << stp1 << "\t" << stp2 << "\t"
1003                                 << avg[0] <<  "\t" << vx_color << "\t" << avg[1]
1004                                 << "\t" << vy_color << "\t" << avg[2] << endl;
1005                     } // end samples
1006                 } // end 2. outer for
1007             } // end 1. outer for
1008             file.close();
1009             qDebug() << "moving" << name << "to safety";
1010             file.copy(name, file_rk_dualbif);
1011             break;
1012         } // end switch
1013     } // end if
1014
1015     if (eul) {
1016         switch(M) {
1017           case 0:
1018           qDebug() << "(euler) standard mode";
1019           qDebug() << "stepsize:" << h;
1020           qDebug() << "number of steps:" << s;
1021           x = t0;
1022           IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1023           k = 0;
1024           hit = 0;
1025           L = 2*PI;
1026           tau = 2*PI/params[4];
1027           while(k < s) {
1028               //qDebug() << k << "<" << s;
1029               //qDebug() << "x is now" << x;
1030               //qDebug() << "deriving values at" << x;
1031               //qDebug() << x0, y0, phi0;
1032               derivs(x,y,dydx,params); // sets first dydx[] for rk4()
1033               derivw(x,y,dW,params,gauss);
1034               euler(y,dydx,dW,x,h,yout);
1035               //for (Int i=0;i<=18;i++)
1036                 //qDebug() << "params["<<i<<"]" << params[i];
1037               if (x >= hit*0.1) {
1038                   //qDebug() << x <<"=="<< hit*0.1;
1039                   //qDebug() << "writing2 to map";
1040                   hit++;
1041                   //qDebug() << "next value to write out" << hit*0.1;
1042                   if (k >= s*fromt) {
1043                       switch(v) {
1044                       case 1:
1045                           euValues_x.insert(x,yout[0]);
1046                       break;
1047                       case 4:
1048                         //xpi = ABS(fmod(yout[0],(2*PI)));
1049                         //ypi = ABS(fmod(yout[1],(2*PI)));
1050                         euValues_x.insert(x,yout[0]);
1051                         euValues_y.insert(x,yout[1]);
1052                         //rkValues_x.insert(x,xpi);
1053                         //rkValues_y.insert(x,ypi);
1054                         euValues_phi.insert(x,yout[2]);
1055                         x1 = yout[0] + params[7]*cos(yout[2]);
1056                         y1 = yout[1] + params[7]*sin(yout[2]);
1057                         x2 = yout[0] - params[8]*cos(yout[2]);
1058                         y2 = yout[1] - params[8]*sin(yout[2]);
1059                         euValues_p1x.insert(x,x1);
1060                         euValues_p1y.insert(x,y1);
1061                         euValues_p2x.insert(x,x2);
1062                         euValues_p2y.insert(x,y2);
1063                         euValues_vx.insert(x,dydx[0]);
1064                         euValues_vy.insert(x,dydx[1]);
1065                         euValues_vphi.insert(x,dydx[2]);
1066                         break;
1067                     } // end fromt switch
1068                   } // end fromt if
1069               } // end hit if
1070               //qDebug() << "making step of size" << h;
1071               for (Int j=0;j<n;j++)
1072                   y[j] = yout[j];
1073               x += h;
1074               //qDebug() << "new x is now" << x;
1075               if (x >= k*tau) {
1076                   //qDebug() << x << ">=" << k*tau;
1077                   //qDebug() << "incrementing k";
1078                   k++;
1079                   //qDebug() << "new value of x needed is" << k*tau;
1080               }
1081           } // end while
1082           qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
1083           qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
1084           qDebug() << "sizeof(euValues_phi) =" << euValues_phi.size();
1085           write_to_file(file_eu_out, euValues_x, euValues_y, euValues_phi, params);
1086           write_to_file(file_eu_pt1, euValues_p1x, euValues_p1y, euValues_phi, params);
1087           write_to_file(file_eu_pt2, euValues_p2x, euValues_p2y, euValues_phi, params);
1088           write_to_file(file_eu_vel, euValues_vx, euValues_vy, euValues_vphi, params);
1089           break;
1090           case 2:
1091           qDebug() << "bifurcation mode";
1092           qDebug() << "bifurcation parameter: " << param1;
1093           qDebug() << "Number of samples : " << samples;
1094           Doub r;
1095           for (Int j=0;(j*stepr1)<=(tor1-fromr1);j++) {
1096               r = fromr1+j*stepr1;
1097               qDebug() << "at step" << j << "of" << (tor1-fromr1)/stepr1;
1098               for (Int z=0;z<samples;z++) {
1099                   x = t0;
1100                   x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
1101                   //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
1102                   IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1103                   bifparam(v,param1,r,params);
1104                   // inner for: iterate over runge-kutta
1105                   k = 0;
1106                   L = 2*PI;
1107                   tau = 2*PI/params[4];
1108                   Doub hnew;
1109                   //qDebug() << "tau: " << tau << "\tL: " << L;
1110                   while (k<s) {
1111                       xtest = x+h; // test next x increment
1112                       //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1113                       if (xtest >= k*tau) {
1114                           //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1115                           //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1116                           hnew = k*tau-x;
1117                           //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1118                           derivs(x,y,dydx,params);
1119                           derivw(x,y, dW, params,gauss);
1120                           euler(y,dydx,dW,x,hnew,yout);
1121                           switch(v) {
1122                           case 1:
1123                               xpi[z][k] = yout[0];
1124                               if (k >= s*fromt) {
1125                                   euValues_x.insert(r,ABS(fmod(xpi[z][k],L)));
1126                                   euValues_xdot.insert(r,dydx[0]);
1127                               }
1128                               break;
1129                           case 4:
1130                               xpi[z][k] = yout[0];
1131                               //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
1132                               ypi[z][k] = yout[1];
1133                               phipi[z][k] = yout[2];
1134                               euDist_x.insert(r,SQR(yout[0]-x0));
1135                               euDist_y.insert(r,SQR(yout[1]-y0));
1136                               euDist_phi.insert(r,SQR(yout[2]-phi0));
1137                               if (k >= s*fromt) {
1138                                   euValues_x.insert(r,ABS(fmod(xpi[z][k],L)));
1139                                   euValues_y.insert(r,ABS(fmod(ypi[z][k],L)));
1140                                   euValues_phi.insert(r,ABS(fmod(phipi[z][k],L)));
1141                                   euValues_xdot.insert(r,dydx[0]);
1142                                   euValues_ydot.insert(r,dydx[1]);
1143                                   euValues_phidot.insert(r,dydx[2]);
1144                               }
1145                               break;
1146                           } // end switch
1147                           k++;
1148                       }
1149                       derivs(x,y,dydx,params);
1150                       derivw(x,y,dW, params,gauss);
1151                       euler(y,dydx,dW,x,h,yout); // xold
1152                       for (Int j=0;j<n;j++)
1153                           y[j] = yout[j]; // y[x+h] !!!!!!!
1154                       x += h; // do it for real
1155                   } // end while
1156                   //qDebug() << "at" << x << "in time for k = " << k;
1157                   switch(v) {
1158                       Int m_x, m_y, m_phi;
1159                       case 1:
1160                       m_x = yout[0]/L;
1161                       avg[0] = (Doub) m_x/k;
1162                       euAvgvelo_x.insert(r,avg[0]);
1163                       break;
1164                       case 4:
1165                       m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
1166                       avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
1167                       trajectory_length = sqrt(SQR(yout[0]-x0)+SQR(yout[1])-y0);
1168                       if (abs(yout[0]-x0) > trajectory_length) {
1169                           velocity_angle = PI;
1170                       } else {
1171                           velocity_angle = (acos((yout[0]-x0)/trajectory_length));
1172                       }
1173                       bias_angle = params[10];
1174                       deflection = (180/PI)*ABS(velocity_angle - bias_angle);
1175                       euAvgvelo_x.insert(r,avg[0]);
1176                       euAvgvelo_y.insert(r,avg[1]);
1177                       euAvgvelo_phi.insert(r,avg[2]);
1178                       euDefAngle.insert(r,deflection);
1179                       break;
1180                   } // end switch
1181               } // samples loop
1182               VecDoub msd_x(samples,0.0);
1183               VecDoub msd_y(samples,0.0);
1184               VecDoub msd_phi(samples,0.0);
1185               Doub sigma_x[s-1];
1186               Doub sigma_y[s-1];
1187               Doub sigma_phi[s-1];
1188               //qDebug() << "samples:" << samples << "steps:" << s;
1189               for(Int i=0;i<s;i++) {
1190                   //qDebug() << "["<<i<<"]:---------------------------------";
1191                   for(Int j=0;j<samples;j++) {
1192                       //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
1193                       msd_x[j] = xpi[j][i];
1194                       //qDebug() << "["<<j<<"]:" << sper_x[j];
1195                       msd_y[j] = ypi[j][i];
1196                       msd_phi[j] = phipi[j][i];
1197                   }
1198                   diffusion(msd_x, samples, sigma);
1199                   sigma_x[i] = sigma;
1200                   diffusion(msd_y, samples, sigma);
1201                   sigma_y[i] = sigma;
1202                   diffusion(msd_phi, samples, sigma);
1203                   sigma_phi[i] = sigma;
1204                   //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
1205               }
1206               // the lists of keys in QMultiMap are read out backwards...
1207               for(Int k=(s-1);k>=0;k--) {
1208                   euVariance_x.insert(r,sigma_x[k]);
1209                   euVariance_y.insert(r,sigma_y[k]);
1210                   euVariance_phi.insert(r,sigma_phi[k]);
1211               }
1212               //qDebug() << euVariance_x;
1213           } // end outer for
1214           /*
1215            * we need to do something about this
1216            */
1217           qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
1218           qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
1219           qDebug() << "sizeof(euValues_xdot) =" << euValues_xdot.size();
1220           qDebug() << "sizeof(euValues_ydot) =" << euValues_ydot.size();
1221           write_to_file(file_eu_bif, euValues_x, euValues_y, euValues_phi, params);
1222           write_to_file(file_eu_biv, euValues_xdot, euValues_ydot, euValues_phidot, params);
1223           write_to_file(file_eu_avg, euAvgvelo_x, euAvgvelo_y, euAvgvelo_phi, params);
1224           write_to_file(file_eu_variance, euVariance_x, euVariance_y, euVariance_phi,params);
1225           write_to_file(file_eu_distance, euDist_x, euDist_y, euDist_phi, params);
1226           write_to_file(file_eu_defangle, euDefAngle, NullMap, NullMap, params);
1227           break;
1228           case 3:
1229           qDebug() << "dual bifurcation mode";
1230           qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
1231           qDebug() << "Number of samples : " << samples;
1232           Doub stp1,stp2;
1233           QString vx_color, vy_color;
1234           QString name = "tmp.dat";
1235           QFile file(name);
1236           QTextStream out(&file);
1237           file.open(QIODevice::WriteOnly);
1238           out << "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl;
1239           for (Int it1=0;(it1*stepr1)<=(tor1-fromr1);it1++) {
1240               for (Int it2=0;(it2*stepr2)<=(tor2-fromr2);it2++) {
1241                   stp1 = fromr1+it1*stepr1;
1242                   stp2 = fromr2+it2*stepr2;
1243                   qDebug() << stp1 << "("<<tor1<<")" << stp2 << "("<<tor2<<")";
1244                   for (Int z=0;z<samples;z++) {
1245                       x = t0;
1246                       x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
1247                       IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1248                       bifparam(v,param1,stp1,params);
1249                       bifparam(v,param2,stp2,params);
1250                       // inner for: iterate over runge-kutta
1251                       Int k = 0;
1252                       L = 2*PI;
1253                       tau = 2*PI/params[4];
1254                       Doub hnew;
1255                       //qDebug() << "tau: " << tau << "\tL: " << L;
1256                       while (k<s) {
1257                           xtest = x+h; // test next x increment
1258                           //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1259                           if (xtest >= k*tau) {
1260                               //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1261                               //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1262                               hnew = k*tau-x;
1263                               //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1264                               derivs(x,y,dydx,params);
1265                               derivw(x,y, dW,params,gauss);
1266                               euler(y,dydx,dW,x,hnew,yout);
1267                               k++;
1268                           }
1269                           derivs(x,y,dydx,params);
1270                           derivw(x,y, dW,params,gauss);
1271                           euler(y,dydx,dW,x,h,yout); // xold
1272                           for (Int j=0;j<n;j++)
1273                               y[j] = yout[j]; // y[x+h] !!!!!!!
1274                           x += h; // do it for real
1275                       } // end while
1276                         Int m_x, m_y, m_phi;
1277                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
1278                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k;
1279                         avg[2] = (Doub) m_phi/k;
1280                         set_color(avg[0], vx_color);
1281                         set_color(avg[1], vy_color);
1282                         out << stp1 << "\t" << stp2 << "\t"
1283                                 << avg[0] <<  "\t" << vx_color << "\t" << avg[1]
1284                                 << "\t" << vy_color << "\t" << avg[2] << endl;
1285
1286                   } // end samples
1287               } // end 2. outer for
1288           } // end 1. outer for
1289           file.close();
1290           qDebug() << "moving" << name << "to safety";
1291           file.copy(name, file_rk_dualbif);
1292           break;
1293         } // end switch
1294     } // end if
1295
1296     if (!ruk && !eul)
1297         throw("No method specified!");
1298
1299     rkValues_x.clear();
1300     rkValues_xdot.clear();
1301     rkValues_y.clear();
1302     rkValues_ydot.clear();
1303     rkValues_phi.clear();
1304     rkValues_phidot.clear();
1305     rkAvgvelo_x.clear();
1306     rkAvgvelo_y.clear();
1307     rkAvgvelo_phi.clear();
1308
1309     euValues_x.clear();
1310     euValues_xdot.clear();
1311     euValues_y.clear();
1312     euValues_ydot.clear();
1313     euValues_phi.clear();
1314     euValues_phidot.clear();
1315     euAvgvelo_x.clear();
1316     euAvgvelo_y.clear();
1317     euAvgvelo_phi.clear();
1318
1319 } // end calculator
1320
1321 //=============================================================================